Naval  Research  Laboratory 

Washington,  DC  20375-5320 


NRL/MR/7140  -05-8835 


Application  of  Differential 
Geometry  to  Acoustics:  Development 
of  a  Generalized  Paraxial  Ray-Trace 
Procedure  from  Geodesic  Deviation 


David  R.  Bergman 

Saint  Peter's  College 
Jersey  City ,  NJ 


January  18,  2005 


Approved  for  public  release;  distribution  is  unlimited. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including 
suggestions  for  reducing  this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis  Highway, 
Suite  1 204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of 
information  if  it  does  not  display  a  currently  valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


1.  REPORT  DATE  (DD-MM-YYYY) 

2.  REPORT  TYPE 

3.  DATES  COVERED  (From  -  To) 

18-01-2005 

Memorandum 

June  2002-August  2004 

4.  TITLE  AND  SUBTITLE 

5a.  CONTRACT  NUMBER 

Application  of  Differential  Geometry  to  Acoustics:  Development  of  a  Generalized 


5b.  GRANT  NUMBER 


Paraxial  Ray-Trace  Procedure  from  Geodesic  Deviation 


6.  AUTHOR(S) 


5c.  PROGRAM  ELEMENT  NUMBER 

61153N 

5d.  PROJECT  NUMBER 


David  R.  Bergman* 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 


5e.  TASK  NUMBER 
UW032-02-43 
5f.  WORK  UNIT  NUMBER 

71-7810 

8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 


Naval  Research  Laboratory,  Code  7140 
4555  Overlook  Avenue,  SW 


Washington,  DC  20375-5320 


NRL/MR/7 140-05-8835 


9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


10.  SPONSOR  /  MONITOR’S  ACRONYM(S) 


Office  of  Naval  Research 
800  North  Quincy  Street 
Arlington,  VA  22217-5660 


ONR 


11.  SPONSOR /MONITOR’S  REPORT 
NUMBER(S) 


12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 


Approved  for  public  release;  distribution  is  unlimited. 


13.  SUPPLEMENTARY  NOTES 

*Saint  Peter’s  College,  Jersey  City,  NJ  07306 

14.  ABSTRACT 

In  this  report,  the  application  of  abstract  differential  geometry  to  geometric  acoustics  is  explored.  The  results  of  this  application  are:  (1)  a 
generalized  paraxial  ray-trace  procedure  valid  for  acoustic  propagation  in  a  random  media  with  subsonic  flow,  (2)  the  demonstration  of  a  con¬ 
tinuum  of  equivalent  paraxial  systems  related  by  a  conformal  transformation,  and  (3)  a  unified  approach  to  treating  problems  in  acoustics,  which 
leads  to  generalized  versions  of  Snell’s  law,  Fermat’s  principle,  and  range-  and  travel-time  integrals  for  layered  media.  The  geodesic  deviation 
vector  is  used  to  model  beam  deformation  and  provides  one  with  an  all-purpose  tool  for  measuring  geometric  transmission  loss  and  locating 
caustics  within  a  ray  skeleton  without  repeatedly  solving  the  ray  equations.  When  applied  to  layered  media,  the  deviation  vector  is  solved  exactly. 
Compared  to  traditional  approaches,  the  results  are  equivalent.  However,  the  difference  is  vast  when  implemented  in  numerical  ray-trace  codes. 
Applications  are  made  to  several  depth-dependent  scenarios,  including  piecewise-linear  sound-speed  and  fluid-velocity  profiles  for  which  the 
exact  caustic  structures  are  determined. 


15.  SUBJECT  TERMS 


Underwater  acoustics;  Propagation;  Geodesic  deviation;  Ray  tracing;  Geometric  acoustics;  Subsonic  flow 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

18.  NUMBER 

19a.  NAME  OF  RESPONSIBLE  PERSON 

OF  ABSTRACT 

OF  PAGES 

Daniel  Wurmser 

a.  REPORT 

b.  ABSTRACT 

c.  THIS  PAGE 

UL 

66 

19b.  TELEPHONE  NUMBER  (include  area 

Unclassified 

Unclassified 

Unclassified 

code) 

(202)  404-4817 

1 


Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  Z39.18 


CONTENTS 


INTRODUCTION . 1 

GENERALIZED  PARAXIAL  RAY  TRACING . 3 

CONFORMALLY-EQUIVALENT  RAY  SYSTEMS . 14 

EXAMPLES . 19 

DISCUSSION  AND  CONCLUSION . 36 

ACKNOWLEDGEMENTS . 37 

REFERENCES . 37 

APPENDIX  A:  CONCEPTS  FROM  DIFFERENTIAL  GEOMETRY . 40 

APPENDIX  B:  NULL  GEODESICS  AND  ACOUSTIC  RAYS . 45 

APPENDIX  C:  CONSTRUCTION  OF  THE  INTERNAL  FRAME  FIELD . 53 

APPENDIX  D:  RIEMANN  TENSOR  AND  CHRISTOFFEL  SYMBOLS 

FOR  THE  ACOUSTIC  METRIC . 57 

APPENDIX  E:  FERMAT’S  PRINCIPLE  AND  SNELL’S  LAW . 59 


iii 


FIGURES 


1 .  Space-time  geodesic  deviation . 5 

2.  Equal-time  deviation . 5 

3.  Beam  deformation . 6 

4.  Initial  wavefront . 11 

5.  Caustic  detail . 20 

6.  Ray  fan  for  canonical  Munk  profile . 20 

7.  Caustic  structure  for  a  harmonic  waveguide . 23 

8.  Convergence/divergence  zones  for  a  quadratic  fluid  duct . 25 

9.  Sample  ray  fan  and  caustic  detail  for  a  quadratic  fluid  duct . 26 

10.  Convergence/divergence  zones  for  linear  sound-speed  and  fluid  profiles . 28 

11.  PlECEWISE-LINEAR  ENVIRONMENTAL  PROFILES .  29 

12.  JACOBI  FIELD,  PIECEWISE  PROFILE .  30 

1 3 .  Ray  fan  incident  on  a  linear  half  space . 34 

14.  Sound-speed  profile . 36 

15.  Caustic  detail  for  a  wedge  duct . 36 

16.  Jacobi  field . 37 

Bl.  Sample  of  a  complete  space-time  ray  structure . 51 

B2.  Relationship  between  fluid  velocity,  ray  tangent,  and  wavefront  normal . 51 


TABLES 

1.  Ray- theory  and  differential-geometry  analogues . 4 


V 


NOTATION 


F3 

M 

xk 

D 

xa 

TP(M ) 
K(M) 
Va 

Va 

{//; «/? } 
D 


U  •  V 

fix  v 
/ 


/' 

V 

(V) 

(Vf 

dA 

d 

A®B 

sJk 


Euclidean  space 
Abstract  manifold 
Cartesian  coordinates  in  Ei 

Vector  in  E3 

Coordinates  in  an  open  patch  of  M 
Tangent  space  at  a  point  p  eM 

Cotangent  space  at  a  point  p  eM 

Components  of  an  element  of  Tp(M),  contravariant  vector 
Components  of  an  element  of  T*(M ) ,  covariant  vector 

Christoffel  symbols  of  the  second  kind 
Christoffel  symbols  of  the  first  kind 
Riemann  curvature  tensor 

Hyperplane  spanned  by  the  vectors  {ai } 

Inner  product  of  vectors  a  and  b  (abstract  notation) 

Dot  product  of  two  vectors  in  £3 
Cross  product  of  two  vectors  in  E 3 

Total  derivative  with  respect  to  the  independent  parameter, 
df  /dt ,  df  Id/ 1 ,  etc.,  depending  on  context 

Spatial  derivative  when / depends  on  only  one  space  variable,  df  jdz 

Gradient  operator  in  E 3 

Column  vector  with  components  8 1 

Row  vector,  (5  x  8  8 ,) 

Covariant  derivative  operator  in  M 

Total  directional  derivative  along  A,  8t  +  A-V 

Total  derivative  in  E J ,  8t  +  r  ■  V 

Direct  product  tensor  of  two  vectors  with  components  AjBl 
Fluid  shear-compression  tensor  8  -uk  +8ku  - 
Fluid  vorticity  8 jVk  ~8kvj 


Fatin  indices,  k,  run  from  1  to  3,  and  Greek  indices,  a  ,  run  from  0  to  3,  with  “0” 
reserved  for  the  time  coordinate.  The  Einstein  summation  convention  is  used  with  sum 
implied  when  any  index  appears  only  twice  as  a  subscript  in  E3,  and  once  as  a  subscript 
and  once  as  a  super  script  in  M,  i.e. 

k= 1  a= 0 


VI 


APPLICATION  OF  DIFFERENTIAL  GEOMETRY  TO  ACOUSTICS: 
DEVELOPMENT  OF  A  GENERALIZED  PARAXIAL  RAY-TRACE 
PROCEDURE  FROM  GEODESIC  DEVIATION 


1.  INTRODUCTION 

When  background  fluid  motion  is  subsonic  the  acoustic  rays  in  a  time-dependent 
moving  ideal  fluid  are  equivalent  to  the  null  geodesics  of  a  pseudo-Riemannian 
manifold  ’  .  This  connection  suggests  a  novel  approach  to  dealing  with  linear  acoustics. 
Equating  the  features  of  ray  theory  with  concepts  from  differential  geometry  offers  an 
advantage  to  traditional  views  by: 

•  Providing  a  unified  approach  to  problems  in  acoustics  that  encapsulates  all  known 
results  in  an  easily  accessible  fonnat. 

•  Leading  to  immediate  generalizations  of  known  problems  that  include  the  effects 
of  random  media  on  the  acoustic  field. 

Specifically,  by  using  geometric  principles  this  approach  leads  to  generalized 
versions  of  Fermat’s  principle,  ray  integrals  for  layered  moving  media,  Snell’s  law,  and 
the  laws  of  geometric  spreading.  In  addition  to  reproducing  known  results,  the 
application  of  differential  geometry  leads  to  a  generalized  version  of  paraxial  ray 
tracing3,4  derived  from  the  law  governing  geodesic  deviation5"9.  The  differential 
geometric  approach  to  ray  tracing  also  offers  a  new  way  to  look  at  different  presentations 
of  ray  theory  as  being  related  by  a  type  of  geometric  invariance  called  conformal 
invariance10.  Conformal  invariance  allows  one  to  view  traditional  approaches  to  ray 
theory  as  belonging  to  a  continuum  of  equivalent  representations. 

Paraxial  ray-trace  procedures  calculate  ray  paths  within  a  ray  bundle  by  solving  a 
system  of  second-order  linear  differential  equations  along  one  ray  path  rather  than 
repeatedly  solving  the  ray  equation  with  different  initial  conditions3,4,11.  From  the  results 
of  this  procedure  the  neighboring  rays  close  to  a  specific  solution  of  the  ray  equations 
may  be  mapped  thus  allowing  one  to  model  the  deformation  of  a  ray  bundle  as  it 
propagates  through  a  random  environment.  This  type  of  ray  tracing  has  become  widely 
used  in  seismology  where  motion  of  the  medium  can  be  ignored. 

The  approach  presented  here  finds  its  natural  home  in  the  time  domain  treating 
acoustic  rays  as  trajectories,  or  histories,  in  space  and  time  in  contrast  to  the  traditional 
view  which  treats  the  rays  as  unit-speed  curves  in  space.  In  the  conventional  approach  a 
ray  trace  is  viewed  as  static  structure  in  space,  sometimes  called  a  skeleton  upon  which 
the  field  flesh  lives.  In  contrast,  the  space-time  approach  views  the  rays  as  encoding 
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information  about  the  history  and  future  evolution  of  the  acoustic  field  as  well  as  its 
shape  at  any  given  instant.  This  infonnation  is  contained  in  the  geometry  of  the  “acoustic 
metric”  and  can  be  thought  of  as  a  property  of  the  medium  itself  for  any  source  placement 
and  geometry. 

Rays  emerge  as  a  high-frequency  approximation  to  the  wave  theory  of  light  or  sound. 
The  theory  of  bicharacteristics  is  more  general,  describing  the  propagation  of  a  surface  of 
initial  data  through  space-time22'23.  The  geodesic  structure  that  detennines  the  rays  in  the 
high-frequency  limit  is  identical  in  form  to  that  of  the  bicharacteristics  of  the  field 
equations. 

The  rate  of  separation  of  neighboring  geodesics  from  a  specified  geodesic  with  similar 
initial  conditions  is  detennined  by  the  stationary  values  of  the  variation  of  the  geodesic 
path  with  respect  to  initial  conditions.  This  variation  is  measured  by  the  Jacobi  field  of 
the  corresponding  geodesic  equation  determined  by  the  second  variation  of  the  action. 
The  same  concept  gives  rise  to  the  notion  of  geodesic  deviation,  a  technique  used  in 
general  relativity  to  describe  gravitational  tidal  forces5'7. 

Most  of  the  applications  presented  in  this  report  make  use  of  the  Jacobi  field  to 
directly  calculate  or  estimate  relevant  aspects  of  an  acoustic  ray  skeleton  and  the  acoustic 
field.  Once  the  Jacobi  field  is  known  for  a  given  situation,  either  analytically  or 
numerically,  it  can  be  used  to  calculate  the: 

1.  Field  intensity  in  the  limit  of  geometric  acoustics. 

2.  Coordinates  of  rays  near  a  pre-chosen  central  ray  within  the  ray  skeleton. 

3.  Exact  location  of  caustics  along  the  central  ray. 

The  behavior  of  the  Jacobi  field  is  detennined  by  the  Riemann  curvature  tensor 
derived  from  the  acoustic  metric — see  Appendix  A  for  details.  Briefly,  a  positive 
curvature  causes  convergence  of  neighboring  geodesics  with  a  common  point  source  and 
predicts  the  existence  of  conjugate  points,  equivalent  to  caustic  formation  in  the  acoustic 
field.  Negative  curvature  leads  to  divergence  of  geodesics  with  a  common  point  source, 
at  roughly  an  exponential  rate,  while  zero  curvature  leads  to  linear  divergence.  For  the 
special  case  of  constant  curvature  the  manifolds  described  by  the  previous  cases 
correspond  to  a  sphere,  pseudosphere,  and  Euclidean  space,  respectively.  The  relation 
between  the  environment  and  the  Riemann  curvature  tensor  allows  one  to  detennine,  by 
inspection  in  many  cases,  whether  or  not  convergence  is  eminent,  and  even  estimate  the 
frequency  of  caustic  formation  without  explicitly  solving  for  the  Jacobi  field.  For 
example,  stationary  fluids  with  depth-dependent  sound  speeds  satisfying  c"(z )  >  0 ,  such 
as  the  canonical  Munk  profile  used  to  describe  the  deep  underwater  sound  channel,  are 
known  to  produce  natural  sound  ducts  giving  rise  to  nontrivial  caustic  structure  in  the 
acoustic  field.  The  curvature  in  such  cases  is  proportional  to  c"(z),  corresponding  to 
locally  spherical  geometry.  The  table  below  provides  a  list  of  terms  from  ray  theory  and 
the  corresponding,  or  relevant,  differential  geometric  quantities.  Each  ray  quantity  is 
either  identical  to,  or  is  derived  from,  the  corresponding  geometric  quantity. 
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Ray  Theory 

Differential  Geometry 

Ray 

Geodesic 

Derivative  of  ray  coordinates 
with  respect  to  initial  conditions 

Jacobi  field.  Geodesic  deviation 

Rays  and  wavefronts  for  a  point  source 

Geodesic  polar  mapping 

Caustics 

Conjugate  points 

Energy  conservation, 

Fermat’s  principle 

Null  constraint 

Ray  integrals  for  layered  media, 

Snell’s  law 

Isometry, 

Symmetry  of  the  metric  tensor 

Focusing  effects  of  the  medium 

Riemann  curvature  tensor 

Table  1.  Ray-theory  and  differential-geometry  analogues 

The  main  purpose  of  this  report  is  to  present  a  complete  set  of  paraxial  ray  equations 
for  three-  and  four-dimensional  acoustic  ray  tracing  in  an  arbitrary  enviromnent, 
including  fluid  motion  and  time  dependence  '  ,  along  with  illustrative  examples 
employing  the  techniques  of  differential  geometry,  and  to  discuss  the  conformally- 
equivalent  representations  of  ray  theory. 

In  Section  2  the  complete  paraxial  ray-tracing  procedure  is  presented  in  four¬ 
dimensional  space-time  and  three-dimensional  space,  along  with  detailed  descriptions  of 
the  quantities  involved.  The  special  case  of  layered  media  leads  to  an  exact  solution  for 
the  Jacobi  field.  Section  3  describes  the  effects  of  conformal  invariance  in  stationary 
fluid  systems  and  compares  the  results  for  four  different  popular  choices  of  ray 
parameter.  Section  4  presents  an  analysis  of  several  systems  that  have  exact  solutions  for 
the  rays  and  the  Jacobi  field,  and  describes  the  caustic  structure  as  predicted  by  the  Jacobi 
field.  Appendix  A  provides  a  list  of  relevant  concepts  from  differential  geometry  as  a 
quick  reference  for  the  reader.  The  reader  interested  in  a  deeper  presentation  should 
consult  Refs.  7,  8,  and  24.  Appendix  B  presents  a  complete  derivation  of  the  geodesic 
structure  associated  with  the  bicharacteristics  beginning  from  first  principles,  comparing 
the  general  case  to  that  of  linear  acoustics  in  the  high-frequency  limit.  Appendix  C 
presents  the  technical  details  behind  constructing  a  ray-centered  basis  for  the  paraxial 
ray-trace  procedure.  Appendix  D  lists  the  components  of  the  Christoffel  symbols  of  the 
second  kind  and  the  Riemann  curvature  tensor  in  a  Cartesian  coordinate  basis,  while 
Appendix  E  introduces  isometry  and  derives  the  ray  integrals  along  with  Snell’s  law  for 
layered  media  with  linear  and  polar  symmetries.  While  Appendices  A,  B,  and  E  are 
mostly  for  completeness  and  rigor,  Appendices  C  and  D  are  an  integral  part  of  the 
paraxial  ray-trace  procedure. 

2,  GENERALIZED  PARAXIAL  RAY  TRACING 

In  this  section  the  complete  generalized  paraxial  ray-tracing  procedure  is  presented  in 
affine  and  time  parameterizations  along  with  results  for  layered  media.  Figure  1 
illustrates  the  basic  features  of  a  ray  system  in  a  two-dimensional  slice  of  space-time. 
The  fiducial  ray  is  labeled  yF  and  neighboring  rays  are  labeled  yN  .  The  coordinates  of 
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each  ray  are  x F  and  x%  respectively,  and  are  related  by  x%  =  xF  +  YM  .  The  choice  of  ray 
parameter  used  in  the  geometric  approach  corresponds  to  an  affine  parameter,  labeled  X  . 
The  significance  of  X  is  discussed  in  Appendix  A. 4.  The  ray  tangent,  7V' ,  in  space-time 
and  the  deviation  vector,  Y“ ,  at  two  values  of  the  affine  parameter,  X ,  are  shown.  The 
deviation  vector  lies  in  a  2-dimensional  space-like  hyperplane  at  each  point  of  the  ray  yF . 
The  Jacobi  equation  describes  the  evolution  of  this  vector  projected  into  this  hyperplane 
and  requires  the  introduction  of  a  ray-centered  orthonormal  basis  (e,\,  I  =  1,  2  (not 

shown  in  the  figure).  In  general,  the  deviation  vector  will  turn  into  the  time  direction 
while  always  remaining  space-like. 


Figure  1.  A  fiducial  geodesic  and  one  Figure  2.  Similar  to  Fig.  1  but  including  the  equal- 

of  its  neighbors  are  shown,  along  with  time  deviation  vector,  shown  at  both  points  along  the 

the  fiducial  tangent  vector  and  the  4-D  fiducial  ray.  The  4-D  and  3-D  deviation  vectors  are 

deviation  vector  at  two  values  of  X ,  initially  equal  by  choice.  As  the  system  is  tracked 

labeled  0  and  1 .  along  yF  ,  Y  picks  up  a  time-like  component  while  Y 

contains  only  space-like  components  at  all  times. 


Coordinate  time  is  a  more  natural  choice  of  ray  parameterization  from  the  physical 
point  of  view.  In  the  coordinate-time-parameterized  version  of  the  problem,  the  deviation 
vector  remains  tangent  to  the  wavefront  in  space  at  all  points  along  the  ray — see  Fig.  2. 
In  this  case,  the  deviation  vector  approximates  a  small  arc  in  the  eikonal  surface  and  may 
be  used  to  reconstruct  the  wavefront  in  a  small  region  about  the  central  ray  by  sweeping 
the  deviation  vector  about  the  central  ray  for  a  complete  set  of  initial  launch  angles, 
forming  a  small  disk  (Fig.  3). 
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Figure  3.  A  3-D  ray  path  is  shown  along  with  a  small  ring  formed  by  the  locus  of  initial  3-D  deviation 
vectors.  As  the  deviation  vectors  are  propagated  along  the  ray  this  ring  deforms  as  a  result  of  the  local 
refractive  index  and  fluid  velocity.  A  pair  of  axes,  labeled  a  and  b,  are  included  to  demonstrate  the  change 
of  orientation  of  the  wavefront. 


2.1  Affine  Representation 

The  paraxial  equations  consist  of  three  main  ingredients:  the  geodesic  equation,  the 
parallel  transport  equation,  and  the  Jacobi  equation,  respectively  listed  below: 


d1 '  xM  „  dxa  dx,s 

■  +  T\p— — —  =  0. 


d/ T 


dX  dX 


(1) 


de?  ,  dxa  p 
~dX  T  aP~dX61  =0, 


d~Y, 


dX 


2  +K1JYJ=  0, 


(2) 

(3) 


where 


y,  =  efg^ 


K  =  R  xax^efiev 

^  IJ  ~  1Vjuccv/3a'  A  • 


(4) 

(5) 


The  explicit  forms  of  YMap  and  Ruavfi  appear  in  Appendix  D.  Many  of  the  longer 
equations  will  be  found  in  Appendices  B-D  along  with  derivations. 
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The  paraxial  ray-tracking  procedure  may  be  outlined  as  follows: 

1.  Find  a  specific  solution  to  the  geodesic  equation,  Eq.  (1),  labeled  the  fiducial 
geodesic,  yF  {A) . 

2.  Solve  the  parallel  transport  equation,  Eq.  (2),  for  the  ray-centered  basis  e“  along 
YfW- 

3.  Using  the  results  of  steps  1  and  2  solve  the  Jacobi  equation,  Eq.  (3),  for  some 
specified  initial  conditions  {or}  to  get  the  components  of  the  deviation  vector  Yr 
along  yF(A)  projected  onto  the  internal  basis. 

4.  The  space-time  components  of  the  deviation  vector,  Yft=Ylef,  give  the 

coordinates  of  a  neighboring  ray,  xF  =  xF  +  . 

Expanding  upon  the  above  procedure  provides  the  necessary  information  for  calculating 
the  geometric  spread  of  the  acoustic  beam  in  space-time. 

5.  Repeat  step  3  by  either: 

a.  Choosing  a  new  set  of  initial  conditions  {<?}  on  T7 ,  such  that 
Yj  (a,  A)  •  Yj  (a,  A)  =  0 . 

b.  Rotating  the  initial  deviation  Y(a)  in  the  (e7)  plane  by  a  small  change  of 
initial  conditions. 

6.  The  infinitesimal  cross-sectional  area  of  the  ray  bundle  in  the  {e,}  plane  is  then 
respectively  given  by: 

a.  f(a)xf(a), 

b.  Y(a)xdY(a)/ 2  . 

The  order  of  the  geodesic  and  deviation  equations  is  reduced  by  introducing  new 
variables  pa  =  xa  and  Pj  =  Yj  and  defining  a  20-dimensional  dynamic  vector, 


u  =  xa  ®pa  ©<  ®Yj  ®PI  ,  7=1,2. 


(6) 


The  complete  paraxial  equations  are  then 


du 

dA 


F(u,u) 


(7) 
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with 


f(u,u)  =  Pa  ©  -r aMrpfJpv  © -r aPrPMe;  ©  pt  ©  -kuYj .  (8) 

The  20  equations  are  subject  to  8  constraints.  Provided  the  initial  data  satisfy  these 
constraints,  the  solution  will  always  satisfy  them  in  principle.  There  are  7  constraints  on 
the  8  components  of  e “  (see  Appendix  C)  and  one  constraint  on  pa .  Implementing  all 
of  the  constraints  reduces  the  dynamical  system  to  12  degrees  of  freedom: 

u  =xk  ®pk  ®t  ®R(a)®  Yj  ®Pj,  (9) 

where  t  is  coordinate  time.  The  corresponding  derivative  vector  is  of  the  form 

F(u,u)  =  pk  ®  Ak  ®w®  a>®  PI  ©  -KuYj ,  (10) 

where  w  =  dt/dA  is  detennined  by  the  null  constraint,  Appendix  B.4,  Eq.  (B.15), 
Ak  =  -Tk oji p" pr' ,  with  all  occurrences  of  p°  replaced  using  the  null  constraint  and  that 
R(a)  is  an  SO(2)  rotation.  Imposing  the  constraints  on  the  system  does  not  necessarily 
lead  to  a  simpler  system.  The  reduction  of  variables  leads  to  an  increased  complexity  in 
the  remaining  equations,  and  the  quantities  co  and  w  typically  contain  denominators  that 
may  vanish  on  occasion. 

An  important  consequence  of  this  construction  is  that  the  deviation  vector  points  in 
four-dimensional  space-time  from  the  fiducial  geodesic,  yF(A),  to  a  point  on  a 
neighboring  geodesic,  yN  (A) ,  with  the  same  value  of  affine  parameter  A  .  This  means 

that  the  deviation  vector  will  not  necessarily  remain  “in  phase”  with  the  fiducial  ray  path 
in  the  traditional  sense  of  the  tenn.  This  does  not  pose  any  problems  in  the  geometric 
description  of  neighboring  curves  since  the  Jacobi  field  simply  maps  out  the  local  space- 
time  geometry  in  a  tube  surrounding  yF(A)  without  prejudice  to  any  coordinate,  time 
being  just  another  coordinate  in  the  four-dimensional  geometric  paradigm.  Two 
adjustments  can  be  made  which  render  the  system  of  Eqs.  (l)-(5)  useful  for 
reconstructing  the  wavefronts  and  predicting  geometric  loss.  These  two  adjustments  are 
independent  of  each  other  and  the  resulting  procedure  is  identical  to  that  outlined  above. 

2.2  Time  Parameterization 

Including  the  auxiliary  basis  described  in  Appendix  C  to  the  above  procedure  allows 
for  the  identification  of  neighboring  points  on  the  wavefront  at  a  given  value  of  A  by 
using  the  data  for  p*1  and  ef  to  construct  eI  (see  Appendix  C.2,  Eq.  (C.3),  for  an 
explicit  form  of  eI ).  Switching  to  a  coordinate-time  parameterization  of  the  problem 
allows  e ,  to  be  tracked  directly  without  the  need  for  e7,  eliminates  the  need  for 
calculating  time,  and  allows  one  to  easily  construct  wavefronts  from  ray  data  by 
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specifying  that  the  ray  equations  be  integrated  to  a  given  set  of  time  values.  The 
corresponding  equations  and  dynamical  vector  are  listed  below: 


*k  =  ^Xk2c2^  Sij  (*<  -  t*j  -°j)+  fa  -  °k  id  0  In  c  +  2r  ■  V  In  c) 


•  Stvk  -cdkc-  cokj  (xj  -  Vj )  ■ +  \  fa.  -  o)kj  ]oj , 


—  e, ,  +  —s!tt(v  xu  +  s  x  77  +  2 Vc  x  n).e, ,  =  0 , 
dt  7’*  2  Jk,y  ,J  Ll 


£L.-kOL+kuYi=  0. 

dt 2  dt 


with  +r0oo  +2r°m0r  and  4  =  (v+V"f"+2V") 


(11) 

(12) 

(13) 


Derivations  of  each  formula  are  given  in  the  appendices.  The  time-parameterized 
representation  of  the  equations  requires  a  16-dimensional  dynamical  array 


u  =  x  ®  p  ®  e)  ®  Yj  ®  Pj , 


(14) 


with  derivative  vector  given  by 

F  =  p®A®  Cie  j®Pj  ©  -KjjYj  .  (15) 

Imposing  constraints  in  the  auxiliary  basis  leads  to  an  1 1 -dimensional  system.  This  is 
not  a  huge  improvement  to  the  12-dimensional  array  of  the  last  section.  However  there 
are  a  few  attractive  features  of  this  representation  of  the  problem.  First,  there  is  no 
equation  for  determining  travel  time  along  the  ray  path  since  time  is  the  parameter  used 
as  a  step  in  the  integration  of  the  system.  Second,  the  equations  are  free  of  any  divisions 
that  might  go  to  zero  during  the  procedure.  Finally,  there  is  the  conceptual  or 
pedagogical  advantage  of  being  able  to  immediately  identify  the  wavefronts  from  the  data 
without  any  additional  work.  The  equations  do  involve  a  greater  number  of  tenns  as 
compared  to  the  affine  counterpart,  introducing  greater  complexity.  The  vector  Yl 
measures  the  deviation  within  the  wavefront. 

2.3  Two  Spatial  Dimensions 

Many  systems  of  interest  in  acoustic  modeling  are  either  legitimately  two-dimensional 
or  can  be  reduced  to  effective  two-dimensional  systems  that  provide  good 
approximations  to  the  true  system.  In  such  cases  the  auxiliary  basis  vector  is  completely 
detennined  at  all  points  along  the  ray  by  a  90-degree  rotation  of  the  wavefront  normal. 
The  reduction  of  spatial  dimensions  also  reduces  the  number  of  components  of  the 
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deviation  vector  from  2  to  1.  In  the  affine  parameterization  there  are  8  degrees  of 
freedom: 


u  =  xa  ®pa  ®Y®P.  (16) 

Moving  to  a  time  parameterization  further  reduces  the  system  to  6  dimensions: 

u  =  x®  p®Y  ®  P .  (17) 

Not  only  is  the  work  required  to  produce  the  answer  significantly  reduced,  but  also  the 
equations  are  usually  simpler  in  form  for  lower-dimensional  systems.  The  basic  steps  in 
the  procedure  are  identical  to  those  outlined  for  the  affine  case. 

2.4  Initial  Conditions,  Mapping  of  Neighboring  Rays,  and  Beam  Deformation 

The  components  of  the  equal-time  deviation  vector  Y,  approximate  a  small  element  of 
arc  length  within  the  eikonal  surface  in  the  eI  direction.  The  initial  value  of  the 
deviation  may  be  chosen  arbitrarily.  The  initial  value  of  Y,  may  also  be  chosen 
arbitrarily  and  is  related  to  the  values  of  the  physical  parameters  that  describe  the 
environment.  In  modeling  the  evolution  of  a  wavefront  the  initial  geometry  is  assumed 
known,  hence  it  is  desirable  to  express  Y10  in  tenns  of  known  quantities.  By  definition 

x^=xF+YM,  YM=x^-Xp,  and  dxF  / dA  =  dtF  / dA{cFnF  +  vF) ,  with  a  similar 
expression  for  the  neighboring  ray,  cF  =  c(xF  )  being  a  shorthand  notation  for  the  local 
sound  speed  evaluated  at  a  point  on  the  fiduciary  ray  (with  similar  expressions  for  all 
other  quantities  evaluated  along  the  ray).  At  the  initial  wavefront  tF  |  =tN  |  by  choice, 

and  the  initial  rates  may  chosen  such  that  dtF  /dA |  =  dtN/dA\Q.  Denoting  this  common 
initial  rate  of  the  time  coordinate  by  the  constant  ft ,  the  initial  velocity  of  the  deviation 

vector  becomes  70  =  ft(cNnN  - cFnF  +0N  -0F),  all  quantities  being  evaluated  at  the 
same  A  but  along  different  rays.  In  tenns  of  time  parameterization 
dY /  dt\  =  cNnN  -  cFnF  +  0 N  -  0F  ,  and  the  initial  deviation  velocity  for  two  rays  with  a 

common  initial  position  dY  /  dt\  =  c0An,  where  A h  =  hN  -hF  with  both  terms  being 

evaluated  at  the  same  point  in  space.  The  quantity  An  determines  the  initial  shape  of  the 
wavefront,  whereas  all  other  quantities  are  assumed  to  be  given  (see  Fig.  4). 
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Figure  4.  The  shape  of  a  small  patch  of  the  initial  wavefront  and  two 
neighboring  wavefront  normal  vectors  are  depicted  for  the  three  cases, 
a:  (a)  concave,  (b)  convex,  and  (c)  planar  section  of  the  wavefront. 


The  initial  values  of  Yt  may  be  expressed  in  terms  of  Y0M .  By  explicitly 
differentiating  the  expression  Yr  =  YMelu  the  following  expressions  for  the  initial  values 
of  Y,  are  derived: 


dY1 
dX  q 


+  -Y'a ) 
2 


ki 


0 


(18) 


where  all  values  are  taken  at  the  initial  point  along  yF(X)  ,  and  dYk  t  dX |  is  given  in  the 
previous  paragraph  in  terms  of  data  at  the  initial  point,  along  both  yF  (X)  and  yN  ( X ) . 
Dividing  by  /?  in  Eq.  (18)  gives  the  appropriate  initial  values  for  dYr  /  dt .  Given  a 
specific  choice  of  initial  wavefront  geometry  at  a  specific  position  in  space-time,  Eq.  (18) 
gives  the  appropriate  initial  conditions  for  dYr  /  dX  or  dY/  /  dt . 


Intuitively  caustics,  or  focal  points,  along  a  ray  may  be  identified  with  points  where 
Y^  =  0  .  To  be  more  explicit  the  deviation  vector  is  expressed  as  a  function  of  the  initial 
conditions  Y(X;(Y0,  Y() ))  =  Yx  (a,) ,  where  ( Y() ,  f(l )  =  a(.  is  a  shorthand  notation  referring  to 
the  set  of  initial  conditions  along  the  7-th  neighboring  ray  to  yF(X)  (with  coordinate 
indices  suppressed).  If  the  initial  conditions  of  the  family  of  neighboring  rays  are  close  to 
those  of  the  fiduciary  ray  (which  is  necessarily  true),  then  a  caustic  point  along  yF  (X)  is 
determined  by  Yx(ai)  =  0.  It  may  happen  that  two  different  neighboring  rays,  say 
yN  (X)  and  yN  ( X ) ,  with  similar  initial  conditions  meet  without  ever  encountering 
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yF  (A)  .  If  the  two  sets  of  initial  conditions  are  close  to  each  other  a  crossing  of  this  pair 
is  determined  by  the  condition  Yk  (a,)  =  YA  ( a2 )  . 


By  applying  step  6b  a  picture  of  the  beam  cross  section  at  any  later  time  may  be 
constructed.  The  area  of  the  planar  cross  section  is  calculated  from 


A(t)  = 


]_ 

2 


Y(t)xdY(t) , 


(19) 


where  Y  and  dY  may  be  expressed  in  terms  of  a.  and  daj.  The  acoustic  intensity  at 
any  point  along  the  ray  is  determined  from  the  cross-sectional  area  of  the  ray  bundle. 
From  the  differential  area  element  dAh  detennined  in  step  5,  one  obtains 

I0/I,  =  ( dAh  ■  t)t  / (dAh  ■  t)0  ,  (20) 


where  t  is  a  unit  vector  tangent  to  the  ray  path  and  Y  is  projected  in  the  wavefront 
normal  basis. 


2.5  Application  to  Layered  Media 

To  find  an  explicit  expression  for  the  deviation  vector  the  tangent  of  the  fiduciary  ray 
and  the  basis  vectors  along  the  ray  must  be  known  at  every  point.  For  two-dimensional 
problems  involving  layered  media  it  is  always  possible  to  find  closed-form  expressions 
for  these  quantities.  To  achieve  this  the  explicit  form  of  the  geodesic  equation,  Eq.  (1),  is 
abandoned  in  favor  of  a  set  of  first-order  equations  derived  from  identifying  isometries  of 
the  metric  tensor  (Appendix  E). 

Consider  a  medium  with  sound  speed  c(z )  and  one-dimensional  fluid  velocity  of  the 
form  ux(z)  =  v(z) .  Rays,  considered  as  curves  in  E' ,  that  are  initially  fired  in  the  x-z 

plane  do  not  turn  out  of  their  initial  osculating  plane,  i.e.  are  torsion  free,  and  constitute 
an  effective  two-dimensional  system.  The  components  of  the  ray  velocity  follow 
immediately  from  the  procedure  in  Appendix  E: 


dt 

dA 


-  va)  , 


(21) 


dx 

dA  ' 


a  +  k, 


0 


(22) 


dz  _  k*0 

dA  c 


- va y 


2  2 

-a  c 


(23) 
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Differentiating  Eq.  (21)  implies  c2t  =  ~{kxv'  +  2cc't)z ,  while  differentiation  of  Eq.  (23) 
and  use  of  the  last  result  implies  z  =  -(/yt/  +  cc'i)i ,  where  dot  denotes  differentiation 
with  respect  to  X  .  It  may  be  demonstrated  that  turning  points  in  time  are  not  allowed,  as 
the  relation  dt  I  d/ 1  =  0  leads  to  the  unphysical  result  {dz  / dX)2  <  0 ,  hence  without  loss  of 
generality,  the  condition  dt  I  d/ 1  >  0  may  be  placed  on  any  ray  in  the  system. 

The  single  basis  vector  for  the  eikonal,  e  ,  is  constructed  by  simply  rotating  the 
nonnal  n  =  dr  /dt  -v  by  ±90  in  the  x-z  plane,  giving  the  following  expressions  for  the 

components  of  the  second-rank  antisymmetric  tensor  ,  exi-e°x  =  ±z/c  and 
ezi-e°z  =  +  kx/c  ,  and  the  useful  relation  exz-e:x  =  ±at0/c  ,  which  is  derived  from  the 
identity  T'ek  -  Tke\  -  Tle,k  -  Tke\  .  Thus  the  explicit  form  of  the  Christoffel  symbols  is 
unnecessary  for  determining  the  ray  tangent  and  its  basis. 

The  Cartesian  components  of  the  covariant  Riemann  tensor  derived  from  the  acoustic 
metric  are  listed  below: 


R 


OxOx 


(24a) 


Ro,x:  =-pr(2cV  +  {V)2  -2  ecu')  , 
4  c 


R, 


OzOz 


=  — ^y(-4c3c"  +  3c2(t>')7  +  4c2uu" -4cc'uur  +  (u'YJ. 


4c 


The  sectional  curvature  along  the  ray  derived  from  the  Riemann  tensor  is 

K  =  -V'(l  - av)-  °^{p')2  -4 c'o'(l  -  au)+—c” , 


which  may  be  reduced  to 


K  =  >^}_d_ 

2  c  dz 


]_d_ 
c  dz 


(a2c2  -(l -ao)2) 


(24b) 

(24c) 


(25) 


The  Jacobi  equation  then  becomes 
d2Y 

-—z-  ±  KY  =  0  .  (26) 

dl2 

The  last  term  in  Eq.  (25)  states  that  c"  governs  the  focusing  properties  of  a  stationary 
medium  with  an  inhomogeneous  sound  speed.  A  ray  propagating  in  a  region  with  c"  >  0 
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will  eventually  encounter  conjugate  points,  while  rays  propagating  in  regions  where 
c"  <  0  will  diverge  from  one  another.  When  the  sound  speed  is  constant  the  effects  of 
fluid  motion  can  be  broken  up  into  two  terms,  aiu"  -(u')2(a/c)2 .  The  second  term  is 
always  negative,  causing  ray  divergence.  In  the  first  term  the  focusing  effects  are 
determined  by  the  relative  sign  of  a  and  v" .  This  term  will  cause  focusing  of  acoustic 
rays  when  a  and  o"  have  the  same  sign.  As  a  consequence  a  constant  fluid  velocity  has 
no  effect  on  the  stability  of  neighboring  rays,  while  a  fluid  velocity  with  a  constant 
gradient  will  interfere  with  the  focusing  caused  by  an  inhomogeneous  sound  speed  with 
c">  0. 

A  more  interesting  situation  occurs  when  both  c  and  u  depend  on  depth.  The 
remarks  of  the  last  paragraph  are  still  valid  with  the  addition  of  an  extra  term  coupling  the 
sound-speed  gradient  to  the  fluid-velocity  gradient,  -t/cV,//c,  appearing  in  Eq.  (25). 
This  tenn  is  more  difficult  to  interpret  than  the  others  contributing  the  curvature. 
Consider  a  simple  situation  in  which  a  waveguide  is  created  by  a  sound-speed  profile 
with  c"  >  0  everywhere.  Above  the  waveguide  axis  c  >  0  ,  while  below  the  waveguide 
axis  c'  <  0  .  If  the  fluid  motion  is  to  the  right  and  characterized  by  a  uniform  gradient, 
then  o'  >  0 .  For  acoustic  rays  fired  to  the  right  a>  0,  resulting  in  a  separation  of 
neighboring  rays  above  the  waveguide  axis  and  an  enhanced  convergence  of  rays  below 
the  waveguide  axis.  When  the  background  fluid  motion  is  weak  and  slowly  varying,  the 
leading-order  terms  in  Eq.  (25)  are  au" ) c2  +  a2  c" )c ,  indicating  that  the  dominant 
effects  are  due  to  the  concavity  of  the  environmental  parameter  values. 

Between  vertical  turning  points  the  deviation  equation  may  be  converted  into  a 
differential  equation  in  the  variable  z  that  in  general  reduces  to  a  single  integral  in  z. 
Changing  variables  to  dcr  =  cdz  leads  to  the  following  general  solution  for  the  Jacobi 
field 


7  =  zcjc1j-p^  +  c2j,  (27) 

with  constants  c\  and  ci  set  to  match  the  initial  conditions  of  the  deviation  vector.  For 
cases  when  z  =  0  identically  along  the  ray,  one  cannot  perform  the  change  of  variables 
from  X  to  z  and  such  cases  must  be  treated  separately.  Equation  (27)  can  be  integrated  in 
some  cases  to  yield  a  closed-form  solution  in  terms  of  ordinary  functions.  For  the  case  of 
layered  stationary  media  Eq.  (27)  becomes 


3 

which  appears  in  a  slightly  different  form  in  Cerveny  . 
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3.  CONFORMALLY-EQUIVALENT  RAY  SYSTEMS 


Although  this  technique  can  be  applied  in  general,  we  now  apply  it  to  two- 
dimensional  underwater  systems  with  range,  depth  and  time  dependence  for  concreteness 
and  simplicity.  The  acoustic  line  element  for  stationary  media  in  four-dimensional  space- 
time  coordinates  is  -  c(x,t)2  dt2  +  dx  ■  dx  =  0.  In  general,  rays  in  space  will  not  be 
confined  to  a  two-dimensional  plane  but  spiral  through  the  medium  as  determined  by  the 
gradients  of  the  sound  speed.  For  situations  in  underwater  acoustics  one  can  find  either 
effective  two-dimensional  subsystems  within  the  ray  structure,  or  in  many  situations  the 
torsion  of  the  rays  is  negligible  for  the  range  of  interest  leading  to 
-  c(z,r,t)2dt 2  +  dz 2  +  dr 2  =  0 ,  with  corresponding  metric  tensor  g uv  =  - c(z,r,t )2  ©  A. . 

As  a  result  of  conformal  invariance  we  are  free  to  pick  and  choose  from  a  continuum 
of  theoretically  equivalent  representations  of  the  same  system.  Why  choose  a  particular 
one  over  another?  Common  to  many  laws  of  physics  is  the  fact  that  the  law  itself  is 
coordinate  or  gauge  invariant  yet  a  specific  representation  must  be  chosen  to  get  the 
problem  solved.  Choosing  a  particular  representation  wisely  results  in  a  simpler  problem 
from  both  the  analytical  and  numerical  points  of  view. 

Below  are  presented  four  popular  choices  of  parameterization  of  the  ray  equations 
along  with  a  comment  of  its  use  in  the  literature.  Each  is  presented  here  as  an  affine 
parameter  leading  to  a  new  geometry.  In  other  references  this  is  taken  to  be  simply  a 
change  in  parameter  devoid  of  any  geometric  significance. 

3,1  Affine  Parameterization 


This  is  the  obvious  choice  from  a  geometric  point  of  view.  It  appeared  once  in  an 
article  by  R.W.  White19,  introducing  the  connection  between  null  geodesics  and  acoustics 
rays.  It  also  appears  to  be  identical  to  the  parameterization  used  by  Cerveny  in  seismic 
ray  theory3.  The  ray  equations  are: 


d2t  1  dt  [  dt  dc 
dA2  c  dA\  dA  dt 


'  dx  dc  dz  dc  3 
yd/ 1  dx  dA  dz , 


=  0, 


d2z 

dc 

(  dt^ 

2 

K 

CN 

^3 

dc 

(  dt^ 

+  c  — 

— 

=  0, 

+  c  — 

— 

dA2 

dz 

1 yd  A) 

dA2 

dx 

\dA  j 

(29) 


(30),  (31) 


The  null  condition  may  be  used  to  eliminate  dt/dA  from  the  equations  for  the  spatial 
variables  if  desired.  Also,  note  that  the  second-order  equation  for  time  may  be  derived  by 
differentiating  the  null  constraint.  When  the  sound  speed  is  explicitly  time  dependent  one 
cannot  perform  the  single  integration  of  dt  =  ds  /  c  unless  the  sound  speed  is  a  separable 
function,  c(t)c(f) .  The  Jacobi  equation  for  this  geometry  is  given  below  in  Cartesian 
coordinates: 
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7  =  0 


(32) 


d2Y 

dA2 


ij 

f  dz  3 

2  ^2 
d  c  i 

f  dx  2 

+  -( 

— 

— 

C1 

yd  A  J 

ax2  1 

ydA  ; 

^  dz  dx  d°  c 
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The  space-time  components  of  the  ray  tangent  are  TM  =  {c  2  c  which  may  look  a 
little  strange  at  first.  This  is  due  to  the  nature  of,  and  units  of,  the  affine  parameter.  One 
can  check  by  hand  that  mn  =  0  as  expected.  The  ray  velocity  infi3  is  given  by 

f  IT0  =  cn  as  expected.  In  ray-centered  coordinates  the  Jacobi  equation  becomes 


d2Y  1  d2c 

^4  +  4^7  =  0 

dA  c 3  dry 


(33) 


Conceptually  this  is  very  easy  to  interpret.  The  stability  of  a  ray  is  governed  by  the 
concavity  of  the  sound  speed  in  the  vicinity  of  the  ray.  In  practice  this  can  be  difficult  to 
implement  as  we  typically  do  not  know  a  priori  c(y)  (sound-speed  data  are  usually  given 
as  c(x,z)).  A  virtue  of  this  parameterization  is  the  simplicity  of  the  equations.  The 
medium  was  not  assumed  to  be  time  independent,  yet  no  time  derivatives  of  the  sound 
speed  appear  in  the  equations  for  either  the  space  coordinates  of  the  ray  or  the  geometric 
spreading. 

3.2  Time  Parameterization 

This  most  closely  resembles  the  equations  found  in  Landau  and  Lifshitz"  .  When  the 
sound  speed  depends  on  only  depth  and  range,  the  time  integral  can  be  performed  along 
the  ray  path  indicating  that  the  differential  equation  is  superfluous.  Time  symmetry  leads 
to  the  constraint 


This  also  indicates  that  travel  time  is  a  good  choice  of  affine  parameter.  By  simply 
rearranging  terms  in  the  line  element  we  can  find  a  representation  of  the  geometry  that 
allows  time  to  be  used  as  an  affine  parameter.  Let  dA  — »  dA/ c2  =  dt  with  k0  =  1 .  The 

new  space-time  metric  structure  is  -  dt2  +  (dx2  +  dz2)/c2  =  0.  The  null  surface  is 
equivalent  to  a  two-dimensional  Riemannian  surface  with  a  line  element  defined  by 
dt 2  =  (dx2  +dz2)/c 2  .  One  can  think  of  this  as  a  travel-time  space  and  the  proper  arc 
length  between  two  points  is  a  measure  of  the  travel  time  along  the  curve  connecting 
them.  In  this  representation  there  is  no  equation  for  time.  The  equations  for  the  ray 
coordinates  are 


2  dxf  dx  dc  dz  dc  2 
c  dt\dt  dx  dt  dz  , 


(34) 
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(35) 


2  dz  f  dx  dc  dz  dc  ' 
c  dtydt  dx  dt  dz  , 


dc 

-c — . 
dz 


The  ray  tangent  in  space-time  is  Tf‘  =  (l  cn)  with  dr/dt  =  cn  .  The  ray  velocity  has 
unit  length  with  respect  to  the  acoustic  metric.  The  Jacobi  equation  for  this  geometry  has 
a  particularly  simple  form 


T  2  y 

— +  jcV2c  -  Vc  •  VcV  =  0  . 

dt~ 


(36) 


Some  attractive  features  of  this  representation  are  that  the  sectional  curvature  does  not 
depend  on  the  ray  momentum  and  time  is  used  as  a  parameter.  Ray  code  based  on  this 
theory  does  not  require  travel-time  calculations  and  wavefronts  can  be  more  easily 
identified  from  numerical  data.  The  ray  equations  are  not  terribly  complicated  as 
compared  with  the  previous  case,  but  the  spreading  equation  is  simpler  in  that  it  does  not 
depend  on  any  of  the  ray-momentum  variables.  The  second  tenn  in  Eq.  (36)  does  lead  to 
theoretically  interesting  changes  and  potential  numerical  effects.  The  curvature  is  a 
defining  property  of  geometry.  As  a  simple  example  consider  the  c-linear  sound-speed 
profile  as  a  function  of  depth.  In  the  affine  parameterization,  this  is  a  space  of  zero 
curvature,  i.e.  flat  as  a  plane.  In  the  time  parameterization,  it  is  a  space  of  constant 
negative  curvature — see  Section  4  for  details.  Additionally,  sound  speeds  with  positive 
concavity  can  have  regions  of  negative  curvature  according  to  Eq.  (36).  In  cases  where 
the  ray  bundles  of  interest  are  required  to  be  parabolic  this  does  not  pose  any  problems. 
In  general,  it  could  lead  to  divergent  behavior  of  Y. 

3.3  Range-Independent  Sound  Speed 

This  choice  of  parameter  is  used  by  researchers  studying  ray  chaos  in  range-dependent 
systems  where  the  ray  equations  are  viewed  as  a  two-dimensional  non-autonomous 
Hamiltonian  system"  ’  .  It  is  important  to  note  that  in  those  cases  a  conformal 
transformation  will  not  make  r  an  affine  parameter.  This  type  of  system  has  also  been 
studied  in  the  context  of  instability  due  to  time-dependent  fluctuations  where  this 
identification  is  possible. 

The  same  trick  used  in  the  last  section  may  be  employed  when  the  sound  speed  is  of 
the  form  c(z,t).  In  this  case,  beginning  from  the  general  case  we  have  one  conserved 
quantity,  dxjdX  =  kx  .  This  indicates  that  range  is  already  an  affine  parameter. 
Rearranging  the  null-line  element  gives  an  effective  two-dimensional  Lorentzian 
geometry,  -dx1  =  - c2dt 2  +  dz2 .  In  this  representation  of  the  ray  structure,  range  is  the 
proper  arc  length  of  the  acoustic  line  element.  Notice  that  the  line  element  is  negative 
indicating  that  the  acoustics  rays  in  ( t ,  z)  space  are  identical  in  structure  to  geodesic  paths 
corresponding  to  massive  particle  trajectories. 
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The  ray  equations  are 


dt  1  dc  (  dt 
— —  - — 

dx  c  dt\dx  j 


7  2  dc  dt  dz 


H - —  0 , 

c  dz  dx  dx 


(37) 


d2z  +  dc  f  dt 
dx2  dz\dxj 


=  0. 


(38) 


The  corresponding  Jacobi  equation  is 


Tr+i^o. 

dx  c  dz 


(39) 


In  this  case  the  curvature  term  is  independent  of  the  ray  momentum.  This  follows  from 
the  reduction  of  the  problem  to  an  effectively  lower-dimensional  system,  and  the  fact  that 
two-dimensional  manifolds  only  have  one  independent  component  of  the  Riemann 
curvature  tensor.  The  equations  take  their  simplest  form  in  this  representation. 

3.4  Three-Dimensional  Arc  Length  Parameterization 

This  is  by  far  the  most  common  presentation  of  the  ray  equations  found  in  texts  on 
acoustics.  It  is  used  in  the  theoretical  treatment  of  ray  geometry  in  E}  as  a  unit-speed 
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curve  and  is  the  basis  for  a  number  of  early  ray  trace  codes  '  . 

In  general,  the  space-time  line  element  implies  the  relation  cdt  =  dS ,  where  S  is 
defined  as  three-dimensional  arc  length  of  the  ray.  This  relation  also  immediately  hints  at 
an  appropriate  confonnal  transfonnation  leading  to  a  geometric  representation  in  which  S 
is  an  affine  parameter.  The  general  form  of  the  complete  set  of  equations  is: 


dt 

~dS 


1 

c 


(40) 


d2x  _  1  dx  f  dc  dz  dc  dx  1  Sc)  1  dc 
dS 2  c  dS  v  oz  dS  dx  dS  c  dt )  c  dx 
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(44) 


One  can  plainly  see  how  complicated  things  get  in  general.  A  significant  feature  of  the 
affine  parameterized  equations,  and  the  1  ©  1  -dimensional  range-parameterized 
equations,  is  that  the  curvature  tensor  only  depends  on  time  through  c  and  not  its 
derivatives.  In  contrast  the  arc-length-parameterized  representation  is  quite  a  bit  more 
complex,  even  for  simple  systems.  Virtues  of  this  representation  are  that  the  ray 
momentum  is  (cos  0  sin  0)  and  the  arc  length  is  known  each  step  of  the  way. 

Sample  ray  traces  using  all  four  approaches  are  given  in  Figs.  5  and  6,  two  examples 
each  for  the  canonical  Munk  profile  and  the  c-quadratic  profile.  A  glance  at  the  detail  of 
a  focal  point  in  the  c-quadratic  case,  Fig.  5,  illustrates  the  precision  with  which  each 
method  located  caustics.  The  different  polygonal  shapes  due  to  different  step-size 
choices  are  pronounced  at  the  turning  points.  Both  sound  speeds  were  expressed  in 
dimensionless  variables  and  the  parameters  in  the  Munk  profile  were  changed  for 
simplicity  to  s  =  0.01 ,  d  =  1.3km  and  d  =  1.5km/s  . 

An  embedded  Runge-Kutta  procedure  with  extrapolation  and  step-size 
control/prediction  was  used  to  integrate  the  system  in  all  four  cases.  This  choice  was 
made  to  take  advantage  of  the  explicit  truncation  error  provided  by  the  output  of  the 
embedded  procedure.  Error  control  was  achieved  by  placing  a  constraint  of  10  6  on  all 
variables  in  the  array,  including  momentum.  This  corresponds  to  an  error  of 
roughly  -  1  jus  in  time  and  - 1  mm  in  distance.  The  null  condition  was  evaluated  at  each 
step  to  gauge  the  relative  success  of  regulating  truncation  error,  rather  than  imposing  a 
tolerance  on  the  constraints  of  the  system.  In  all  cases  the  magnitude  of  the  error  in  light 
cone  placement  ds  =  0  was  bound  from  above  by  -  10  5 ,  and  in  most  cases  by  ~10^6 . 
Caustics,  surface  and  bottom  hits,  as  well  as  potential  detector  interactions  were 
identified  by  a  change  in  sign  of  an  appropriate  quantity,  i.e.  caustics  were  located  by  a 
change  in  sign  of  Y.  The  location  was  obtained  by  the  method  of  false  position. 
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Figure  5.  Detail  of  caustic  for  C  =  1  +  0.0  lcf  :  (a)  range  parameterized  and  (b)  arc-length  parameterized. 


Figure  6.  Sample  ray  trace  for  Munk  profde:  (a)  affine  parameterized  and  (b)  time  parameterized. 


4.  EXAMPLES 

In  this  section  several  examples  having  exact  solutions  to  the  ray  equations  and  the 
Jacobi  equation  in  terms  of  ordinary  functions  are  presented  and  discussed.  Many  of  the 

13172131 

cases  presented  here  may  be  found  in  the  literature  ’  ’  ’  and  have  been  included  here 
to  shed  light  on  the  techniques  outlined  in  Section  2.  The  approach  presented  here  makes 
explicit  use  of  the  paraxial  equations,  Eqs.  (3)  and  (28),  with  attention  paid  to  the 
prediction  of  caustic  curves  from  zeros  of  the  Jacobi  field,  and  interpreting  the  stability  of 
neighboring  rays  by  mapping  divergence  and  convergence  zones  based  on  the  sectional 
curvature. 
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In  principle  one  knows  7(A) .  The  caustic  is  the  locus  of  points  determined  by  7  =  0. 
Solving  this  equation  gives  Ar(a) ,  the  value  of  the  ray  parameter  at  the  caustic  location 
in  terms  of  the  ray  initial  conditions.  Using  the  ray  coordinates  then  gives  the  caustic 
curve  in  space  xc(a),  zc{a),  parameterized  by  the  initial  conditions  of  the  ray.  Cusps 

are  determined  by  finding  values  of  a  for  which  the  velocity  of  the  caustic  curve 
vanishes: 


dxc 

da 


da 


(45) 


A  few  simple  cases  are  analyzed  from  the  differential  geometric  point  of  view.  The 
simplest  cases  are  well  known  examples  that  provide  more  insight  into  the  language  of 
differential  geometry.  Afterwards  less  trivial  cases  are  presented,  with  solutions  to  the 
caustic  curves  determined  by  the  geodesic  deviation  equation. 

4.1  Layered  Media  with  Smooth  c(z)  and  u(z) 


The  first  two  examples  illustrate  the  differences  in  using  A,  and  time  as  affine 
parameters.  Both  cases  are  described  by  stationary  media,  and  the  conformal  metric 
introduced  earlier  is  used.  Recall  from  Section  3.2  that  for  purely  two-dimensional 
systems,  the  sectional  curvature  is  K  =  cV2c-Vc-Vc.  For  simplicity  consider  rays  in 
the  x-z  plane.  Along  any  ray  an  orthonomial  basis  may  be  constructed  by  inspection 
consisting  of  the  unit  tangent  vector  along  the  ray  t  =  xi  +  zk  and  two  basis  vectors,  the 
first  being  obtained  by  rotation  of  the  ray  tangent  by  90  degrees,  ex  =  -zi  +  xk  ,  the  other 
being  e2  =  cj  .  In  the  conformal  metric,  coordinate  time  is  a  proper  affine  parameter  and 
“dot”  refers  to  differentiation  with  respect  to  time.  The  basis  vectors  ( eve2 )  are 
orthonormal  with  respect  to  the  conformal  inner  product  g..  =  Sjj/c2.  The  sectional 
curvature  in  the  (t ,ex)  plane  is  K  =  cc"-(c')2,  while  the  sectional  curvature  in  the 
(f,e2)  plane  vanishes.  From  these  results  it  follows  that  rays  fired  from  a  point  source  at 
the  same  initial  angle  but  in  different  planes  will  diverge  linearly  in  time. 

Spherical  Geometry  and  Ideal  Focusing 

By  direct  calculation  the  sound-speed  profile  c(z)  =  Ccosh(^A^z/c)  produces  a 
space  of  constant  positive  curvature  K  =  K0.  The  deviation  vector  can  be  solved  exactly 

for  all  rays  in  this  system:  Y (t)  =  A  cos(J~Kt)  +  BAn(J~Kt) .  One  can  conclude  from  this 
that  ideal  periodic  focusing  will  occur  along  each  and  every  ray  with  the  same  frequency 
regardless  of  the  placement  of  the  source.  The  time-versus-depth  integral  can  be 
calculated  in  this  case  leading  to  the  result  that  the  travel  time  between  turning  points  is 
independent  of  initial  conditions,  hence  the  point-like  focusing  demonstrated  in  this  case 
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will  always  occur31.  In  the  affine  representation  the  sectional  curvature  is  K  =  a2K0. 

The  appearance  of  the  ray  parameter  leads  to  a  nontrivial  change  in  the  interpretation  of 
the  geometry  of  space  and  its  effect  on  geodesics. 

With  time  parameterization,  the  geometry  is  identical  to  a  two-dimensional  sphere. 
The  periodic  ray  structure  can  be  mapped  into  a  sphere  by  morphing  the  rays  in  such  a 
way  that  they  resemble  those  of  a  source  placed  on  the  waveguide  axis,  then  identifying 
consecutive  focal  points  with  the  north  and  south  poles  of  the  sphere.  The  rest  of  the 
space  is  then  projected  stereographically  onto  the  sphere.  The  resulting  space,  while 
being  geometrically  identical  to  a  sphere,  is  topologically  equivalent  to  a  sphere  with  two 
holes  pinched  in  it  identified  with  z  =  ±oo  . 

With  affine  parameterization,  the  geometry  is  spherical  along  each  ray  with  a  different 
effective  radius  of  curvature  for  each  value  of  a  .  Rays  fired  towards  z  =  ±oo  never  turn 
and  their  neighbors  drift  away  from  them  at  a  constant  rate.  In  the  previous  case  all  rays 
experience  the  same  curvature,  indicating  that  even  those  fired  towards  z  =  ±oo  converge 
with  their  neighbors  passing  through  caustics.  It  can  be  easily  verified  by  direct 
calculation  of  the  ray  integrals  that  the  total  travel  time  from  -  oo  to  +  oo  is  equal  to  the 
period  between  consecutive  conjugate  points.  Hence  for  a  source  placed  anywhere  in  the 
space,  the  vertical  rays  meet  the  end  of  their  journeys  before  ever  encountering  their 
neighbors.  The  finite  travel  time  from  -oo  to  +oo  is  an  indication  that  the  space 
described  in  the  first  case  suffers  from  geodesic  incompleteness  . 

The  Pseudosphere 

The  exact  ray  paths  for  the  sound  speed  c(z)  =  a +  bz  are  circular  arcs.  The  induced 
metric  space  is  known  in  the  literature  as  a  Lobaschevsky  space  and  is  an  example  of  a 
space  of  constant  negative  curvature,  or  pseudosphere,  K  =  -b 2 .  The  deviation  vector  in 
this  case  may  be  found  by  inspection:  Y(t)  =  Acos\\{bt)  +  f?sinh(bt) ,  implying  that  the 
distance  between  neighboring  rays  grows  exponentially  in  time.  Once  again  the  affine 
version  of  the  theory  predicts  a  slightly  different  result.  The  curvature  in  this  case  is 
identically  zero  indicating  that  rays  diverge  linearly  in  X  . 

Harmonic  Waveguide 

The  harmonic  waveguide,  described  by  n2  =a  +  bz2 ,  has  been  the  subject  of  many 
studies  due  to  the  existence  of  exact  solutions  with  nontrivial  behavior  for  the  rays  and 

modes  associated  with  this  profile '°.  Consider  the  special  case  c_1  =  z)2  ,  with 

the  units  of  speed  scaled  to  1 .  This  sound  speed  forms  a  waveguide  with  symmetry  axis 
at  z  =  1  and  c  — »  oo  at  z  =  0  and  z  =  2.  The  exact  ray  paths  are  sinusoidal: 

I -  g 

z  =  l  +  -v/l  —  p2  sin( A  +  cpn ) ,  x  =  pA,  sin  <p0  =  , 0  (46) 

VI  -p2 
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Consider  a  point  source  placed  on  the  waveguide  axis,  z  =  1 .  The  sectional  curvature  is 
positive  definite,  the  corresponding  Jacobi  field  is 


_  /?2sinA  +  (l-/f)AcosA 
■yj\  -  (1  -  p2)sm2  A 


(47) 


For  the  special  cases  p  =  1  and  p  =  0,  7  =  sin  A  and  7  =  A,  respectively,  as  expected. 
The  caustic  is  determined  by  solving  the  transcendental  equation 

tan  A  = - £-A.  (48) 

P 

A  plot  of  the  rays  and  the  caustics  for  a  source  on  the  symmetry  axis  is  given  in  Fig.  7. 
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Figure  7.  Ray  fan  and  caustic  structure  for  the  n  quadratic  sound-speed  profile.  The  rays  are  plotted  for 
even  increments  of  cos<90  ■  The  caustics  are  traced  by  inserting  solutions  to  Eq.  (48)  into  Eq.  (46). 


The  next  few  cases  involve  fluid  motion.  Very  few  interesting  cases  with  exact  solutions 
exist.  We  consider  first  the  curvature  induced  by  a  class  of  symmetric  velocity  profiles 
that  are  found  to  create  ducted  regions.  These  ducts  are  identified  by  mapping  the 
regions  of  positive  sectional  curvature  into  the  x-z  plane.  Afterwards,  the  exact  ray  paths 
are  presented  for  two  special  cases  with  nontrivial  sound  speed  and  fluid  velocity. 
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Wind-Induced  Sound  Ducts 


In  general,  effects  of  o"  on  a  ray  system  are  illustrated  for  a  class  of  velocity 
functions  of  the  form  u(z)  =  c0{z/L )2" ,  where  c0  is  the  local  sound  speed,  assumed  to  be 
a  constant  L  ,  a  length  scale  chosen  such  that  the  flow  becomes  supersonic  for  |z|  >  L  and 

n  >  0  is  a  positive  integer.  This  class  of  functions  has  the  following  common  properties: 
u(z)  is  a  monotonically  increasing  function  for  z  >  0,  u(z)  >  0  for  all  z,  v '(z)  >0  for  z  > 
0,  o"(z)  >  0  for  all  z,  and  o(-z)  =  v(z) .  For  the  time  being  attention  will  be  paid  to  rays 
that  remain  in  regions  of  subsonic  flow. 

From  Snell’s  Law,  £(z)  +  sec<9  =  s0  +  sec#0,  where  s0  and  #0  are  the  Mach  number 
and  initial  angle  of  inclination  at  the  source  point,  respectively.  For  purposes  of 
illustration  it  is  assumed  that  the  source  is  place  at  z  =  0  and  attention  focused  on  rays 
fired  in  the  upper  half  plane.  For  rays  fired  against  the  wind,  #0  >  njl  and  sec#0  <  -1 . 
Combined  with  Snell’s  law  this  information  implies  that  sec#  <  0  at  all  points  along  the 
ray,  indicating  that  rays  with  a  negative  ray  parameter  will  never  turn  back  toward  the 
source.  On  the  other  hand  rays  with  positive  ray  parameter  will  always  encounter  vertical 
turning  points,  and  because  of  the  symmetry  of  u(z)  this  will  happen  periodically. 


The  vertical  turning  points  may  be  found  from  Snell’s  law  by  setting  #  =  0  and 
solving  for  z.  Although  the  form  of  Snell’s  law  used  here  is  expressed  in  terms  of  the 
direction  of  the  wavefront  normal,  the  geometry  of  the  rays  indicates  that  0  =  0  when  the 
ray  encounters  a  vertical  turning  point.  The  equation  for  the  turning  points  becomes 
s(z)  =  -l  +  sec#0.  The  condition  2>sec#0  =^>  0O  <n/2  is  imposed  to  single  out  rays 
that  turn  before  entering  supersonic  regions,  a  condition  that  holds  for  any  choice  of  o . 
The  sectional  curvature  for  this  class  of  problems  is 


K  = 


k: 


2n{4n-\)C2n-2  cos#, 


t-2  2 

L  c0 


COS  #n 


2  n 


cos  #  An  -  1 


(49) 


for  all  rays  in  the  system,  where  the  dimensionless  variable  C,  =  z  / L  is  used.  Clearly  the 
sign  of  K  depends  on  two  factors,  cos#0  and  the  term  in  parenthesis.  Equation  (49)  will 
change  sign  when  cos#  =  cos#0(2 - l/2n) ,  the  same  result  holding  for  rays  fired  both 
into  and  against  the  wind.  Since  cos#  is  an  increasing  function  along  rays  with 
cos  #0  <  0 ,  the  condition  is  never  satisfied  and  these  rays  will  remain  in  a  divergence 
zone  forever.  Rays  with  cos  #0  >  0  are  a  little  more  interesting.  The  curvature  will 
remain  positive  along  these  rays  only  if  cos#  <  cos#0(2  - 1  / 2/z) .  Since  v  is  strictly 
increasing  cos  #  is  strictly  increasing  as  well  and  it  suffices  to  check  the  curvature  at  the 
turning  point  where  cos  #  =  1 .  The  curvature  will  be  zero  at  the  turning  point  if 
cos#0  =  2n/(4n-V) .  For  n  =  1  this  leads  to  the  constraint  cos#0  =2/3,  or  #0  «48°. 
Rays  fired  into  the  wind  with  cos  #0  >2/3  will  remain  in  a  region  of  positive  curvature 
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as  they  propagate  downstream.  Rays  with  cos#0<2/3  will  start  with  K  >  0  and 
eventually  along  the  ray  K  <  0 ,  causing  neighboring  rays  to  begin  diverging  with  respect 
to  A,  in  general  when  cos#  =  2ncos#0  /(4n  - 1) .  Figure  8  shows  the  convergence  and 

divergence  zones  for  a  point  source  placed  on  the  waveguide  axis.  The  zones  are  labeled 
by  the  sign  of  the  sectional  curvature,  K. 

In  the  case  n  =  1,  the  ray  equations  and  the  deviation  equation  may  all  be  integrated 
between  turning  points  to  yield  solutions  for  the  time  of  flight,  range,  and  deviation 
vector  as  functions  of  depth  in  terms  of  elliptic  functions.  A  ray  trace  using  the  same 
technique  employed  in  Section  3  is  shown  below  along  for  one  half  of  a  ray  fan,  with  the 
placement  of  the  caustics  along  each  ray  shown  in  the  detail  (Fig.  9).  The  procedure  was 
stable  with  an  upper  bound  on  the  deviations  from  the  null  surface  ~10  4  for  the  same 
truncation-error  tolerance.  Figure  8  shows  the  regions  of  positive  and  negative  curvature 
mapped  in  the  x-z  plane  for  a  source  at  the  origin.  There  is  a  cut-off  for  rays  fired  with 
the  wind  at  p  =  2/3,  z  =  1/V2  .  Rays  with  p  >  2/3  have  turning  points  at  z  <  1/  and 
remain  completely  within  a  region  in  positive  curvature. 


Figure  8.  Convergence  and  divergence  zones  for  the  quadratic  fluid-velocity  profile.  Regions  are 
marked  according  to  the  sign  of  K.  The  outer  K  =  0  boundary  corresponds  to  the  <2  =  0.  The  inner 
boundaries  correspond  to  a  source  placed  at  the  origin,  on  the  symmetry  axis  of  the  guide. 
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Figure  9.  Ray  trace  for  (left)  the  quadratic  wind  profde,  with  detail  (right)  for  one  cycle  of  the  ray- 
path  motion  illustrating  caustic  formation. 


Linear  Fluid  Velocity 

For  a  linear  fluid-velocity  profile,  v  =  cz/ L,  c  =  constant,  and  K  =  - k j2 / L1  .  The 
deviation  equation  may  be  integrated  to  yield 

Y  =  {yJ4k  )sinh(  VZ/t)  +  70  cosh(VZ/l) .  (50) 

Hence  neighboring  rays  with  the  same  initial  position  will  diverge  at  a  rate  of 
Y  =  (to / vZ)sinh(«/y  2 / c)  with  respect  to  the  affine  parameter.  The  ray  coordinates 
may  also  be  solved  in  this  case.  Defining  new  dimensionless  variables  along  the  ray 
£  =  a(c £)  -1,  £  =  z/L,  r  =  c2at/L ,  ri  =  xlL,  <j  =  ah/L  ,  and  the  new  ray  parameter 
(5  =  ca  ,  Eqs.  (21)-(23)  may  be  integrated  to  yield 


C,  =  P  cosh(cr  -  <p) 

(51) 

r  =  P  sinh  (cr  -  (p)  -  A0 

(52) 

£  =  £o  +^lnA  +^cr--^(2-^cosh(cr-^))sinh(cr-^)  , 

(53) 

where  the  parameters  A0  =  -  fi1  +  <^0  and  (p  =  \n(j3/A0)  have  been  defined.  Eqs. 

(5 1)-(53)  represent  the  affine -parameterized  coordinates  for  a  segment  of  the  ray  with 
z  >  0  and  the  initial  conditions  set  so  that  r(0)  =  0 .  The  ray-path  geometry  implied  by 

1  T 

these  equations  is  well  known  : 
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Equation  (54)  gives  the  upward  part  of  the  ray  while  the  downward  portion  is  determined 
by  symmetry  to  be  2 tj  -  rf ;  tj0,  g(l  are  the  initial  coordinates  of  the  ray  and  rj  ,  gw  are 

the  coordinates  of  the  turning  point. 

Linear  Fluid-Velocity  and  Sound-Speed  Profiles 

Generalizing  the  last  case,  we  consider  the  case  v  =  £  and  c  =  1  +  .  The  exact  ray 

paths  are 


A  =  - 


a2(s2- 1) 


1(1  -  avf  -  a  c  + 


o2(s2  -l)3 


-arctan 


(1  -  av)  +  cas 
/(I  -  av)2  -  a2c2  Va 


The  sectional  curvature  is  K  =  - a  (a  +  s)/c3 .  In  this  case  the  curvature  is  negative  for 
rays  fired  with  the  wind,  a  >  0  ,  indicating  divergent  behavior.  For  rays  fired  against  the 
wind  the  sign  on  the  curvature  depends  on  the  sign  of  (a  +  s) .  Setting  (a  +  s)  =  0 
implies  cos6*0  =  —s  .  For  s  >  1  the  fluid  motion  is  always  subsonic  and  a  change  in  sign 
of  K  for  rays  fired  against  the  wind  is  not  possible,  while  for  1  >  s  >  0  fluid  flow 
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becomes  supersonic  at  z  =  1/(1  -  s) ,  K  >  0  when  |cos#0|  <  s  and  K  <  0  when  Icos  <90 1  >  s  . 
The  curve  defined  by  a  =  0  is 


7a=0 


1 


e  U  +  ^o  j 


(57) 


and  the  curve  defined  by  a--s  when  1  >  e  >  0  is 


V \~£2  -  ,  1 

Va=-s  = - (#“4)- 


-ln 


^  l  +  £^ 


s24\-£2  U  + 


(58) 


A  plot  of  the  convergence  zones  for  each  case  is  given  in  Fig.  10. 


*1  n 


Figure  10.  Convergence  and  divergence  zones  for  the  c-linear,  V -linear  case  outlined  above. 

(a)  corresponds  to  £  =  0.5  and  is  sound-speed  dominated  with  subsonic  flow,  while  in  case 

(b)  £  =  2  ,  with  supersonic  flow  occurring  at  %  =  1 . 


4.2  Piecewise-Linear  Profiles 

Although  the  linear  sound-speed  and  fluid-velocity  profiles  lead  to  ray  systems  with 
somewhat  benign  behavior,  they  are  frequently  used  to  fonnulate  exact  ray  theoretical 
models  of  real  situations  where  the  sound-speed  and  velocity  profiles  are  approximated 
by  piecewise-linear  functions.  In  such  cases  the  ray  paths  and  their  first  derivatives  are 
continuous  and  exact  solutions  to  the  ray  path  and  the  Jacobi  field  are  known  in  each 
piecewise-linear  region.  The  discontinuity  in  the  sound-speed  and  velocity  gradients 
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leads  to  the  presence  of  a  Dirac  delta  function  in  the  sectional  curvature,  in  turn  leading 
to  special  boundary  conditions  for  matching  the  Jacobi  fields  of  different  segments  at  the 
boundary  between  two  regions.  The  basic  situation  is  illustrated  in  Fig.  11.  The  medium 
is  divided  in  depth  into  a  finite  number  of  regions.  As  a  ray  travels  through  the  medium 
it  will  turn  and  periodically  encounter  each  region,  in  principle  an  indefinite  number  of 
times.  The  ray  segments  are  indexed  in  order  of  increasing  ray  parameter  without 
reference  to  the  corresponding  region.  The  Jacobi  field  is  a  function  of  the  affine 
parameter  evaluated  along  the  ray  path.  Each  ray  segment  in  space  corresponds  to  an 
increment  of  the  affine  parameter.  Passing  from  one  region  to  another  in  space 
corresponds  to  boundary  point  along  the  A  axis — see  Fig.  12. 


Range 


Figure  11.  A  physical  situation  approximated  by  piecewise-linear  sound-speed  and 
fluid-velocity  profiles.  A  sample  ray  is  drawn  illustrating  the  meaning  of  the  index 
labels  along  the  ray  as  well  as  those  labeling  the  different  regions. 


Once  the  rays  are  matched  up,  the  parameter  value,  time  of  flight,  and  range  for  each 
segment  may  be  determined.  The  Jacobi  field  is  continuous  leading  to  the  boundary 
condition  Yi+l=Yr  Integrating  Eq.  (26)  across  a  boundary  leads  to  the  following 

condition  on  the  discontinuity  in  Y : 
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YM(AB)-Y1(AB)  +  QYI(KB)  =  (), 


(59) 


Q  =  -  2  (aAt>'(l  -  auB)  +  a2cBAc')  ,  (60) 

^BCB 

where  the  subscript  B  means  evaluated  at  the  boundary  and  A  o',  Ac'  are  the 
discontinuities  in  the  fluid  velocity  and  sound  speed,  respectively,  at  the  boundary.  Once 
these  conditions  are  matched  for  a  significant  number  of  segments  the  caustics  and  their 
features  are  determined  in  the  same  manner  as  before.  A  striking  feature  of  these  models 
is  that  nontrivial  focusing  emerges  as  a  result  of  the  boundary  conditions,  exclusively  the 
curvature  being  otherwise  zero  or  negative.  In  general  the  rays  and  the  Jacobi  field  in 
each  region  are  given  by  Eqs.  (55)  and  (56)  of  the  last  section  with  appropriate  changes  in 
scale  of  the  coordinates. 


Figure  12.  The  Jacobi  field  for  the  profile  illustrated  in  Fig.  11.  The  individual  segments 
of  Y  shown  here  are  linear  merely  for  simplicity.  The  different  depth  regions  in  Fig.  11 
correspond  to  intervals  on  the  A  axis.  Comparing  Figs.  12  and  11  indicates,  for 
example,  that  the  4th  and  6th  intervals  in  Fig.  12  both  correspond  to  region  3  of  Fig.  11. 


Two  special  cases  of  interest  are:  (1)  a  stationary  medium  with  a  piecewise-linear 
sound  speed  and  (2)  constant  sound  speed  with  a  piecewise-linear  fluid  velocity.  In  case 
1  the  curvature  is  identically  zero  in  each  region  and  the  Jacobi  field  may  be  written 
Y:  =  A:  +  BjK ,  while  in  case  2  the  curvature  is  a  negative  constant  in  each  region  and  the 
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Jacobi  field  is  Yt  =  Aje'r~AA  +  B {e  ^A .  In  each  case  applying  boundary  conditions  on  Y 
leads  to  expressions  for  the  coefficients  A  and  B  in  terms  of  the  ray  parameter  p  and  the 
source  coordinates.  The  caustic  in  each  region  is  determined  by  the  condition  Yi  =  0 

leading  to  AiC ,  with  AiC  =  -Ai/Bi  in  case  1  and  2y/-KAiC  =  ln(-  Bj  At)  in  case  2,  for 

the  value  of  affine  parameter  along  the  7-th  segment  of  the  ray  at  which  the  caustic 
occurs.  Inserting  this  expression  into  the  appropriate  ray  equation  then  gives  the  caustic 
coordinates  along  the  ray  as  a  function  of  initial  conditions.  The  piecewise-linear  profiles 
give  rise  to  cusp  and  broken  caustics  depending  on  the  relative  values  of  the  profile 
slopes  in  consecutive  regions.  Two  cases  utilizing  this  procedure  are  illustrated  next. 
Both  cases  involve  stationary  media  due  to  their  simplicity  and  past  use  in  underwater 
acoustic  modeling.  The  generalization  to  non-stationary  media  based  on  the  steps 
outlined  above  is  immediate. 


Sound  Incident  on  a  Linear  Half  Space 


This  classic  problem  of  sound  in  a  homogeneous  medium  incident  on  a  region  of 
space  with  a  graded  refractive  index  is  known  to  produce  a  caustic  in  the  emerging 
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field  ,  and  in  the  graded  medium.  In  affine  parameterization,  the  equations  for  the  ray 
path  are 


dx 

dA 


dt  _  1  dz  _  Vl -a2c2 

dA  c2  ’  dA  c 


(61) 


with  the  local  sound  speed  given  by  c(z)  =  co0(z) +  c0(l  -  z/Z)0(-z) ,  where  0(z)  is 
the  Heaviside-step  function  and  the  dimensions  of  A  are  L7T.  Defining  cos6*0  =  p  the 
rays  in  region  z  >  0  are  straight  lines  given  by 

r/  =  pA  +  r/0,  £  =  4  ±  V1-  P 2  A  >  r  =  A  +  r0 .  (62) 


Here  we  define  the  dimensionless  variables  A  =  A/Lc0,  £  =  z/L ,  rj  =  x/L,  and 
r  =  tc0  /  L  .  For  simplicity  we  set  r0  =  0  and  consider  a  source  placed  at  r/0  =  0  and  £0 
>  0.  Furthermore,  we  focus  attention  on  rays  with  sin  0Q  <  0.  The  solution  of  the  ray 
equations  in  the  linear  region  may  be  expressed  as 

(f7  - (f7i  +  P  1 V1  - P2 ))  +(£“1)2  =  P  2  >  (63) 

where  77,  =  cot^ .  The  equation  for  77  is  the  same  in  both  regions.  The  ray  paths  in 
space  are  circles  of  radius  r  =  |sec6*0|  centered  at  (77,  +tan  0O  l).  By  symmetry  the  ray 
path  emerges  from  region  2  at  -0O.  Hence  the  third  segment  satisfies  Az/Ax  =  -  tan  <90 . 
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At  the  second  boundary  point  we  have  77  2  =  <^0  cot  60  +  2  tan  0O ,  leading  to  the  following 
equation  for  £ : 


«A)  =  1/irprA-(^+2(p-2-l))  .  (64) 

The  initial  conditions  on  Y]  are  chosen  to  describe  the  geometry  of  neighboring  rays 
from  a  point  source,  Yx  =  86  A .  Applying  the  boundary  conditions  leads  to  the  following 
recursion  relation 


^(A)  =  ^_1(A)-/^_1(Ai._1)(A-Am),i  =  2,3  ,  (65) 

for  the  solution  along  each  segment  of  the  ray  path,  with  k  =  p2  / ^1- p2  .  Solving  for 
the  coefficients  gives  A2=  1  -  kAx  ,  A3  =  -1 ,  B2  =  /cA] ,  and  B3  =  4 k1  ,  with 

72  =(l-/tA1)A  +  /tA:j,  Y3=-A  +  4k-\  (66) 

for  the  Jacobi  field  in  each  region.  An  overall  factor  of  86  has  been  divided  out  for 
convenience.  For  a  caustic  to  exist  in  the  lower  portion  of  the  medium  it  is  sufficient  to 
require  F2  ( A  2)  <  0  .  This  places  the  following  constraint  on  the  initial  conditions  of  the 

ray  path:  >  2 (/T2  -l).  This  condition  implies  two  things  depending  on  the  question 

asked.  First,  for  a  given  initial  ray  labeled  p,  at  what  location  of  the  source  does  the 
caustic  on  that  ray  pop  out  (in)  the  lower  half  space?  For  gn  >  2(p  2  -l),  the  caustic 

along  ray  p  will  be  inside  the  lower  region.  Asked  another  way,  for  a  fixed  source 
location  which  rays  pass  through  the  caustic  in  region  1  and  which  in  region  2?  Solving 
for  p  in  the  inequality  gives 


P  > 


(67) 


Rays  fired  at  steep  angles  (small  p)  will  penetrate  deeply  and  pass  through  the  caustic 
in  the  upper  half  space  while  those  with  shallow  launch  angles  (large  p)  will  pass  through 
the  caustic  in  the  lower  half  space.  These  statements  assume  that  the  ray  in  question  will 
pass  through  a  caustic. 

Along  ray  segment  3,  we  have  Ac  =4 k  1 .  Using  the  ray  coordinates  gives 


&  =2(p  2 -l)-£o  ’  0c  =4^p~2-l 

or  simply 


(68) 
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(69) 


for  the  caustic  curve  in  the  upper  half  space.  The  caustic  curve  in  the  lower  portion  of 
space  is  more  difficult  to  extract.  For  a  caustic  to  form  along  segment  2  of  the  ray  path 

we  have  Ac  =  /tA2j/(/tA,  -1) .  To  extract  a  curve  f(^c,Jjc)  =  0  for  the  caustic, 

independent  of  the  initial  conditions,  requires  solving  the  following  equation  for  p: 

[a2  +  b2)p6  -b(2b  +  \)pA  +  (2b  +  \)p2  -1  =  0  ,  (70) 

where  b  =  \  +  g{)  and  a  =  gg  / .  The  exact  roots  may  be  found  by  substituting  p2  =  q 

and  solving  the  third-order  equation.  Rather  than  solve  for  !gc  (r/c )  in  region  2,  we 
express  the  caustic  as  a  parameterized  curve  (gr(  p),/]i:(  p))  for  fixed  gtl : 


-i’ 


(f-i  f=p 


( 

1 

v 


a -p2y 

fe+i  )P2-tf ; 


(71) 


The  velocity  of  the  caustic  is 

djh  ei  pbp  !(£q+3>-3) 

dp  ( } 

_ pip  2('go+3>-  3) _  (73| 

dp  °  fo  +  l)p2-l)s  V+/(^+2fo-2)  +  l-2f0)' 

Clearly  both  vanish  when  p  =  ^3/ (^0  +  3)  indicating  the  development  of  a  cusp.  This 
cusp  will  always  exist  in  region  2  as  long  as  >  0  .  Ray  traces  and  caustics  are  shown 

for  a  few  source  placements  in  Fig.  13.  The  caustics,  appearing  on  top  of  the  ray  trace, 
were  derived  from  the  solutions  presented  here.  One  can  clearly  see  the  development  of 
the  cusp  in  the  lower  half  space.  As  the  source  is  pulled  up,  the  caustic  sinks  further 
down  into  the  lower  region  with  the  cusp  moving  further  away  from  the  origin.  The  tail 
of  the  caustic  approaches  zero  as  rjc  — >  00  .  In  the  limiting  case  where  the  source  in 
placed  on  the  rj  axis,  the  cusp  moves  right  up  to  the  source  and  the  tail  merges  with  the 
77  axis. 
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Figure  13.  Ray  fan  incident  on  a  c-linear  half  space  for  3  source  placements,  from  upper  left  going 
clockwise,  E,  =1,5,  and  10. 


The  Wedge  Duct 

Consider  c  (£)  =  (1  +  £)0(£)  +  (1  -  <j£)&(-£)  (Fig.  14).  The  upper  portion  of  the  duct 
is  fixed  at  a  45-degree  slope,  while  a  measures  the  steepness  of  the  lower  part  of  the 
duct  relative  to  that  of  the  upper  part.  For  simplicity  consider  a  source  placed  on  the 
positive  E,  axis.  The  initial  ray  parameter  is  p  =  cos  O0  /(I  +  ) .  The  ray  is  divided  into 

segments  labeled  yi .  The  points  where  the  ray  crosses  the  //  axis  are  labeled  //,  and  the 
corresponding  parameter  value  A, .  The  first  segment  is  given  by 
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(74) 


(v  +  p-'^l-p2  (!  +  4)2)  +(l  +  %Y=p-2, 

with  a  -  sign  for  rays  fired  upward  and  a  +  sign  for  rays  fired  downward.  After 
traversing  this  segment,  the  ray  periodically  crosses  the  rj  axis  with  the  same  value  of  ray 
parameter  determined  by  Snell’s  law,  cos6(  =  cos#0/(l  +  g0).  In  region  I  the  ray  is 

(r/-(r/i+P^^^-p2]j  +(c^'  =  p2cr^2,  (75) 

while  rays  in  region  II  are  given  by 

(7  -  (7 +  P~l^~P2]j  +  (l  +  #)2  =  P2  •  (76) 

The  crossing  point  of  the  first  segment  is  r/x  =  p  yl  -  p  ±  p  y  1  -  p  (l  +  £0)  ,  and  the 
crossing  points  for  the  other  segments  are 

Vi  =  Vi  +(/ (l  +  a^)-2)p ~lyll-p2  ,  i  =  2,4,  6,... 

Vi  =Vi+  (i  ~  l)(l  +  &  '  )p  X  V1  ~ P2  ’  i  =  3,5,7,... 

for  the  lower  and  upper  arcs,  respectively.  Applying  boundary  conditions  leads  to  the 
following  expressions  for  the  coefficients  and  Bt : 

1  ±  (/  - 1)(1  +  cT1)/?  /  =  1, 3,5,... 

'  _ [- <7(1  ± (i - 1)(1  +  <T~' )R)  i  =  2,4, 6,... 

_  (1^)  ( .  _  i)(i  ±  ( .  _  i)(i  +  ^  1  +  ^ 2 } 

B.=\  aK 

(1  +  ^  ((f  - 1)  ±  (i(z  -  2)  +  2  +  i(i  -  2)<r~1  )/?  +  (i  -1)R2)  1 =  2’4’6”-' 

.  /f 


where,  for  simplicity,  we  have  defined  the  ratio 


7?  = 


A 


and  k  =  p2 /t]{-  p'  as  before. 
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Sample  rays  for  a  =  5  and  £0=1  are  shown  in  Fig.  15  along  with  caustics.  The 
caustics  are  plotted  for  q  e  [0.2, 0.9] .  For  q  =  1  the  upper  and  lower  caustics  terminate  at 
ranges  of  4.15  and  2.08,  respectively.  Both  become  unbound  as  q  approaches  0.  The 
Jacobi  field  along  the  first  three  segments  is  displayed  in  Fig.  15  for  the  four  rays  in  Fig. 
16. 


Figure  14.  Sound-speed  profile  for  the  wedge  duct.  The  lower  portion  is  plotted  for  G  =  5  . 


Figure  15.  Caustic  for  the  bilinear  duct.  Only  4  sample  rays  were  traced  for  cleanliness. 
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Figure  16.  The  Jacobi  field  plotted  for  the  first  3  segments  of  4  different  rays,  <7  =  5  and  g()  =  1  . 


5.  DISCUSSION  AND  CONCLUSION 

The  space-time  dynamic  ray-tracing  procedure  outlined  in  Section  1  provides  a 
generalization  to  the  approach  currently  used  in  seismology,  and  offers  a  method  of  beam 
tracing  in  a  random  medium.  The  identification  of  acoustic  rays  with  null  geodesics  of  a 
pseudo-Riemannian  manifold  leads  immediately  to  this  generalization.  The  geodesic 
deviation  vector,  used  to  model  the  deformation  of  an  acoustic  beam,  allows  one  to 
measure  geometric  transmission  loss  and  caustic  location  without  repeatedly  solving  the 
ray  equation.  Additionally  this  identification  has  led  to  a  continuum  of  theoretically 
equivalent  paraxial  ray-trace  systems  connected  by  conformal  transfonnations,  and  points 
to  a  unified  approach  to  extracting  well-known  results  such  as  Snell’s  law  and  Fermat’s 
principle  in  random  fluid  media.  The  null  hypersurface  is  viewed  as  a  manifold  created 
by  the  local  sound  speed  and  fluid  velocity  which  encodes  all  information  about  the 
acoustic  field. 
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The  development  here  was  intended  for  describing  acoustics  in  the  presence  of 
subsonic  flow.  This  requirement  is  necessary  for  the  identification  of  acoustic  rays  with 
null  geodesics.  For  cases  involving  supersonic  flow,  the  bicharacteristic  condition 
defines  the  geometry  of  the  acoustic  rays  as  Euclidean  line  elements  of  zero  length.  Such 
spaces  are  a  special  case  of  a  wider  class  called  Finsler  spaces.  This  identification 
indicates  that  these  techniques  may  be  generalized  to  include  cases  involving  supersonic 
flow. 

ACKNOWLEDGEMENTS 

The  work  was  performed  under  the  Office  of  Naval  Research  /  American  Society  for 
Engineering  Education  (ONR/ASEE)  Summer  Faculty  Program.  In  addition  to  the  ASEE, 
the  author  would  like  to  thank  the  Naval  Research  Laboratory  (NRL)  in  Washington,  DC 
for  hosting  this  fellowship  for  three  summers  (2002-2004),  during  which  time  some  parts 
of  this  work  were  completed.  The  author  would  also  like  to  acknowledge  Dr.  Daniel 
Wurmser,  NRL,  for  helpful  comments  and  discussions  during  the  start  of  this  project,  and 
Dr.  A.  Lewis  Licht,  University  of  Illinois  at  Chicago,  for  bringing  Matt  Visser’s  article  to 
his  attention.  Finally,  thanks  to  Drs.  Robert  Gragg  and  Roger  Gauss,  NRL,  for  helpful 
discussions. 


REFERENCES 

1.  W.G.  Unruh,  “Experimental  black  hole  evaporation?,”  Phys.  Rev.  Lett.  46,  1351- 
1353  (1981). 

2.  M.  Visser,  “Acoustic  black  holes:  Horizons,  ergospheres,  and  Hawking 
radiation,”  Classical  and  Quantum  Gravity  15,  1767-1791  (1998). 

3.  V.  Cerveny,  Seismic  Ray  Theory  (Cambridge  University  Press,  Cambridge,  2001). 

4.  R.  Guenther,  Modern  Optics  (John  Wiley  &  Sons,  New  York,  1990). 
(Specifically,  chapter  5,  section  entitled  “ Propagation  in  a  Graded  Index  Optical 
Fiber",  for  a  waveguide  paraxial  ray  approximation.) 

5.  A.  Abromowicz  and  W.  Kluzniak,  “Epicyclic  Orbital  Oscillations  in  Newton’s 
and  Einstein’s  Dynamics,”  General  Relativity  and  Gravitation  35,  69-77  (2003). 

6.  S.L.  Bazanski,  “Kinematics  of  relative  motion  of  test  particles  in  general  relativity 
(1),”  Ann.  Inst.  Henri  Poincare,  Section  A,  XXVII,  145-166  (1977). 

7.  C.W.  Misner,  K.S.  Thorne,  and  J.A.  Wheeler,  Gravitation  (Freeman,  New  York, 
1973). 


37 


8.  B.  O’Neill,  Elementary  Differential  Geometry  (Academic  Press,  San  Diego, 
1966),  Ch.  VII. 

9.  P.  Schneider,  J.  Ehlers,  and  E.E.  Falco,  Gravitational  Lenses,  Astronomy  and 
Astrophysics  Library  (Springer-Verlag,  New  York,  1992). 

10.  R.M.  Wald,  General  Relativity  (The  University  of  Chicago  Press,  Chicago,  1984). 

1 1 .  M.  Porter  and  H.  Bucker,  “Gaussian  beam  tracing  for  computing  ocean  acoustic 
fields,”  J.  Acoust.  Soc.  Am.  82,  1349-1359  (1987). 

12.  G.S.  Heller,  “Propagation  of  acoustic  discontinuities  in  an  inhomogeneous 
moving  liquid  medium,”  J.  Acoust.  Soc.  Am.  25,  950-951  (1953). 

13.  E.T.  Kornhauser,  “Ray  theory  for  moving  fluids,”  J.  Acoust.  Soc.  Am.  25,  945- 
949(1953). 

14.  C.  Perkeris,  “Theory  of  propagation  of  sound  in  a  half  -space  of  variable  sound 
speed  under  the  condition  of  formation  of  a  caustic,”  J.  Acoust.  Soc.  Am.  18,  295- 
315  (1946). 

15.  O.V.  Rudenko,  A.K.  Sukhorukova,  and  A.P.  Sukhorukov,  “Full  Solutions  to  the 
Equations  of  Geometrical  Acoustics  in  Stratified  Moving  Media,”  Acoustical 
Physics  43,  339-343  (1997). 

16.  R.J.  Thompson,  “Ray  theory  for  an  inhomogeneous  moving  medium,”  J.  Acoust. 
Soc.  Am.  51,  1675-1682  (1972). 

17. 1.  Tolstoy  and  C.S.  Clay,  Ocean  Acoustics,  Theory  and  Experiment  in  Underwater 
Sound  (McGraw-Hill,  New  York,  1966). 

18.  P.  Ugincius,  “Acoustic-ray  equations  for  a  moving,  inhomogeneous  medium,”  J. 
Acoust.  Soc.  Am.  37,  476-479  (1965). 

19.  P.  Ugincius,  “Ray  acoustics  and  Fermat’s  principle  in  a  moving  inhomogeneous 
medium,”  J.  Acoust.  Soc.  Am.  51,  1759-1763  (1972). 

20.  R.W.  White,  “Acoustic  ray  tracing  in  moving  inhomogeneous  fluids,”  J.  Acoust. 
Soc.  Am.  53,  1700-1704  (1973). 

21.  E.R.  Franchi  and  M.  Jacobson,  “Ray  propagation  in  a  channel  with  depth- variable 
sound  speed  and  current,”  J.  Acoust.  Soc.  Am.  52,  316-331  (1972). 

22.  R.  Courant  and  D.  Hilbert,  Methods  of  Mathematical  Physics:  Volume  II  (John 
Wiley  &  Sons,  New  York,  1962). 


38 


23.  F.  John,  Partial  Differential  Equations  (Courant  Institute  of  Mathematical 
Sciences  Lecture  Notes,  New  York  University,  New  York,  1952). 

24.  S.  Kobayashi  and  K.  Nomizu,  Foundations  of  Differential  Geometry:  Volumes  I 
and  II  (John  Wiley  &  Sons,  New  York,  1963).  ( See  Volume  II,  specifically 
Chapter  VIII for  a  discussion  of  Jacobi  fields  and  conjugate  point  theorems .) 

25.  L.D.  Landau  and  E.M.  Lifshitz,  Fluid  Mechanics,  Second  edition  (Butterworth  & 
Heinemann,  Oxford,  1987). 

26.  K.  Smith,  M.  Brown,  and  F.  Tappert,  “Ray  chaos  in  underwater  acoustics,”  J. 
Acoust.  Soc.  Am.  91,  1939-1949  (1992). 

27.  K.  Smith,  M.  Brown,  and  F.  Tappert,  “Acoustic  ray  chaos  induced  by  mesoscopic 
oceanic  structure,”  J.  Acoust.  Soc.  Am.  91,  1950-1959  (1992). 

28.  H.  Davis  et  al.,  “The  Hudson  Laboratories  Ray  Tracing  Program,”  Technical 
Report  No.  150,  Hudson  Laboratories  of  Columbia  University,  1968. 

29.  L.B.  Palmer  and  D.M.  Fromm,  “The  Range-Dependent  Active  System 
Performance  Prediction  Model  (RASP),”  NRL/FR/5 160-92-9383,  Naval 
Research  Laboratory,  Washington,  DC,  1992. 

30.  T.  Foreman,  “A  Frequency  Dependent  Ray  Theory,”  Technical  Report,  ARL-TR- 
88-17,  Applied  Research  Laboratories,  The  University  of  Texas  at  Austin,  1988. 

31. L.M.  Brekhovskikh,  Waves  in  Layered  Media  (Academic  Press,  New  York, 
1960). 

32.  S.W.  Hawking  and  G.F.R.  Ellis,  The  Large  Scale  Structure  of  Space-Time 
(Cambridge  University  Press,  Cambridge,  1973). 


39 


Appendix  A 


CONCEPTS  FROM  DIFFERENTIAL  GEOMETRY 


A.l  Metric 

The  metric  is  a  generalization  of  the  dot  product  and  is  usually  introduced  via  a 
quadratic  form 

ds2  =  ajjdxldx .  ,  (A.l) 

with  ds  a  differential  arc  length  and  ai/(xk)  continuously  differentiable  functions  of 
coordinates.  One  can  always  find  a  small  neighborhood  about  an  arbitrary  point  such  that 

a„  ~  diag(\..  A  -1...-1)  ,  (A.2) 

n  m 


modulo  constant  scale  factors,  with  n  +  m  =  dim(M). 

When  m  =  0,  the  metric  reduces  to  the  identity  matrix  and  the  geometry  is  locally 
Euclidean.  When  I n  —  m\  =  dim(M)  -2 ,  the  metric  reduces  to  one  of  two  choices  for  the 

Minkowski  metric  used  in  special  relativity,  specifically  for  dim(M)  =  4.  Manifolds  with 
metric  defined  by  Eq.  (A.l)  that  are  locally  Euclidean  are  referred  to  as  Riemannian, 
while  those  that  have  local  Minkowski  geometry  are  termed  either  pseudo-Riemannian  or 
Lorentzian. 

A  consequence  of  Eq.  (A.2)  is  that  curves  on  a  pseudo-Riemannian  manifold  can  have 
positive,  negative,  or  zero  length.  The  coordinates  of  the  manifold  are  divided  into  x, ,  i 

=  1 . .  .dim(M)  -  1,  and  t,  and  choosing  m  =  1.  With  this  convention,  curves  with  ds2  >  0 
are  called  space-like,  and  curves  with  ds2  <  0  time-like.  These  two  types  of  curves  are 
separated  by  a  hypersurface  defined  by  the  class  of  curves  with  ds2  =  0  referred  to  as  a 
null  hypersurface,  a  generalization  of  the  light  cone. 

For  subsonic  flow,  the  metric  introduced  in  Appendix  B  describes  a  Lorentzian 
manifold  with  n  =  3  and  m  =  1 .  The  characteristic  condition  identifies  the  acoustic  rays 
with  geodesic  curves  on  the  null  hypersurface  similar  to  photon  trajectories  in  general 
relativity. 
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A.2  Manifolds 


A  manifold  is  an  abstraction  from  the  theory  of  surfaces  that  does  not  include  a 
surrounding  Euclidean  space  of  higher  dimension  in  its  description.  For  a  detailed 
definition  of  a  manifold,  see  the  References.  We  present  here  only  a  few  concepts  that 
are  used  in  other  places  within  the  text.  A  tangent  space,  T  (M) ,  is  defined  at  each  point 

on  the  manifold  which  serves  as  an  abstraction  of  the  tangent  plane  to  a  surface.  Each 
tangent  space  has  the  structure  of  flat  Euclidean  space,  modulo  constant  scale  factors. 
The  metric  in  each  tangent  space,  g  (xa  )  ,  is  fixed  but  is  allowed  to  vary  from  point  to 

point  in  M.  If  U,V  e  T  (M)  then  their  inner  product  is  given  by  gflvUflVv ,  and 
gMVVf,Vv  =  \vf  defines  the  magnitude  of  a  vector.  The  Vp  are  the  components  of  the 
vector  in  some  chosen  basis  of  T  ( M ) .  A  dual  tangent  space,  T*  (M) ,  is  also  defined  at 
each  point  of  M.  The  metric  tensor  defines  a  mapping  from  Tp  (M)  — >  T*  (M) , 
Tp  =g/JVTv.  The  Tp  are  the  components  of  the  dual  vector  in  some  chosen  basis  of 
T*(M )  .  Tensors  of  arbitrary  rank  are  elements  of  a  direct  product  space  constructed  of 
multiple  copies  of  the  tangent  and  dual  tangent  spaces 
T  (M) x  —  T  (M) x T* (M) x---xT* (M) .  For  example  the  metric  presented  in  Section 

B.l  are  components  of  a  tensor  in  T*  x  T* .  These  spaces  do  not  require  a  concept  of 

distance  to  be  defined.  In  some  sense  defining  a  metric  before  this  section  is  like  putting 
the  cart  before  the  horse.  The  components  of  a  tensor  of  arbitrary  rank  are  denoted 
Ta'"'a* Pv~pM  ■  A  standard  choice  of  basis  for  Tp(M )  and  T*(M )  are  8u  and  dx" 

respectively,  where  xM  are  coordinates  defined  on  an  open  set  of  M  containing  p.  In  this 
basis  a  vector  and  its  dual  are  denoted  V  =  Vfldu  and  VpdxM ,  respectively. 

A.3  Covariant  Differentiation  and  Parallel  Transport 

The  covariant  derivative  of  a  vector  and  dual  vector  are  defined  by 

dmuv  =dpu''+rvMu\  d PUV  =dpuv  - r\vu \  ,  (A.3) 

respectively,  where 

TV  =^gvP(SMgz/3  +  dAgM/}-dfigM).  (A.4) 

These  Christoffel  symbols  take  into  account  the  turning  of  the  basis  vectors  as  one 
differentiates  the  vector,  requiring  comparison  of  V  at  two  points  of  M.  Eq.  (A.3)  is  a 
component  form  of  the  result. 
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Comparing  vectors  at  different  points  in  M  is  not  a  trivial  operation  since  they  live  in 
different  spaces.  To  compare  V  e  T AM)  to  V  e  T AM)  ,  p,q  <=M  with  p  *  q  ,  one  of 
the  vectors  must  be  moved  into  the  space  of  the  other.  In  general  this  requires  specifying 
the  path  in  space,  C,  traversed  by  V  from  p  to  q.  If  T M  is  the  velocity  of  C  then  V  is 
parallel  transported  along  C  if 


T^D/V  =  0  .  (A. 5) 

A.4  Affine  Parameter 

Geodesic  curves  parallel  transport  their  velocity  vector,  i.e. 

DTV 

- -  TMD  Tv  =  0  .  (A. 6) 

d; i  M 

When  Eq.  (A. 6)  holds  the  parameter  X  defined  along  the  curve  is  called  an  affine 
parameter.  Any  other  choice  of  parameterization  will  change  Eq.  (A. 6)  to 

TMDJv=aTv,  (A. 7) 

where  u  is  an  arbitrary  scalar  function  defined  along  the  curve.  Equation  (A. 6)  makes 
the  transported  tangent  vector  “parallel”  to  the  tangent  of  the  curve  at  the  new  point 
(allowing  a  change  in  magnitude),  whereas  the  auto-parallel  condition  requires  the 
transported  tangent  vector  to  be  identical  to  the  tangent  vector  at  the  new  point  along  the 
curve.  For  space-like  geodesics,  arc  length  is  an  affine  parameter.  Time-like  geodesics 
may  also  be  parameterized  by  their  length  interpreted  as  a  standard  internal  time,  or  the 
proper  time  of  an  observer  whose  world  line  coincides  with  the  particular  time-like  curve. 
Clearly,  for  null  curves  arc  length  cannot  be  used  as  a  parameter.  In  such  cases  one 
imposes  Eq.  (A. 6)  and  the  null  constraint  thus  defining  the  parameter  as  affine  with  no 
physical  significance  attached. 

A.5  Geodesics 

Defined  in  the  previous  section  as  auto-parallel,  the  geodesic  equation  is  also  the 
Euler-Lagrange  equation  derived  from  a  variation  of  the  arc  length, 

S\ds  =  0  ,  (A.  8) 

with  ds  defined  in  Appendix  A.l.  Equation  (A. 8)  defines  geodesics  as  curves  of 
“optimal”  length,  either  minimizing  or  maximizing  the  length  between  two  points  of  M. 
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A.6  Curvature 


The  Riemann  curvature  tensor,  and  the  covariant  Riemann  tensor  are  derived  from  the 
Christoffel  symbols: 

R\vp  =  dvT\P  -  dpr\v  +  r^vr\p  -  tmaPt\v  ,  (A.9) 

RmvP  =  dv{p-,ap}-d/}{p-,av}  +  Tpav{p-,pp}-Tpap{p-,pv}.  (A.10) 

The  tensor  defined  in  Eq.  (A.  10)  has  several  symmetries  that  reduce  the  total  number  of 
components: 

n  D  _  _ D  _  _ D 

1Xpvafi  LXa(3pv  LXvpa(5  pv ft  a  * 

The  effects  of  curvature  are  most  easily  seen  by  parallel  transporting  a  vector  around  a 
closed  loop  in  M.  If  this  task  were  attempted  on  a  plane  or  in  Euclidean  space  the  result 
would  be  trivial,  the  transported  vector  would  be  identical  to  the  original  vector.  A 
simple  nontrivial  example  of  parallel  transport  on  a  curved  manifold  is  constructed  by 
moving  a  vector  from  the  north  pole  of  a  sphere  around  an  equilateral  geodesic  triangle 
with  0  =  90° .  The  transported  vector  is  rotated  by  90  relative  to  the  old  vector. 

A.7  The  Jacobi  Equation 

The  geodesic  equation  gives  curves  of  optimal  length  on  M  satisfying  Eq.  (A. 5).  The 
stability  of  a  geodesic  relative  to  slight  variation  in  initial  conditions  is  governed  by  the 
Jacobi  equation.  For  convenience  we  define  the  Jacobi  vector,  Y,  along  the  geodesic 
subject  to  the  constraint  gpvYpxv  =  0 .  Focusing  exclusively  on  the  null  geodesics  of  M 

confines  Y  to  the  space-like  plane  (e, ) .  Geometrically  Y  points  from  one  geodesic  to  a 

point  on  a  neighboring  geodesic  with  the  same  value  of  A  .  The  vector  Y  obeys  a  second- 
order  linear  differential  equation  along  the  geodesic: 

~jjf-+K/jYj  =o,  (A.H) 

which  is  derived  from  the  second  variation  of  the  length  integral.  F7  =  e“ gapYp  is  the 
projection  of  Y  onto  the  geodesic-centered  frame  field  defined  in  Appendix  C,  and  the 
curvature  matrix  Ku  =  Rliavf,xu xfl ep e'j  is  introduced.  For  I  =  J,  Ku ,  no  sum  on  /,  is  the 

sectional  curvature  in  the  plane  spanned  by  (i,  e,  'j ,  and  measures  the  relative  rate  of 
separation  of  neighboring  geodesics. 

Equation  (A.  11)  governs  the  stability  of  the  geodesics  and  provides  a  convenient 
paradigm  for  describing  the  focusing  and  spreading  of  a  geodesic  bundle  in  M.  The 
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tendency  for  neighboring  geodesic  to  converge  or  diverge  is  determined  by  the  matrix 
Ku  .  If  the  curvature  is  positive  definite  in  every  plane  containing  the  ray  tangent,  then 

neighboring  rays  from  a  common  source  point  in  all  directions  tend  to  converge  and  a 
three-dimensional  ray  bundle  will  focus.  If  the  curvature  is  zero  or  negative  everywhere 
along  the  ray,  then  neighboring  rays  from  a  common  source  point  will  tend  to  eventually 
diverge.  (If  Y0  ^  0  then  the  initial  shape  of  the  wavefront  may  cause  the  formation  of  a 
caustic  independently  of  the  focusing  properties  of  the  medium,  as  happens  with  light 
reflected  from  a  concave  surface  in  Euclidean  space.)  The  conditions  Kn  >0,  and 
det  Kjj  >  0  along  y(X)  are  necessary  and  sufficient  for  the  quadratic  form  KuYjYj  to  be 
positive  definite. 

A.8  Conjugate  Points 

Two  points  on  a  geodesic  y(X) ,  P  and  Q,  are  said  to  be  conjugate  if  there  exists  a 
nontrivial  Jacobi  vector  along  y(X)  that  vanishes  at  P  and  Q.  According  to  the 
conjugate-point  theorems  of  differential  geometry,  if  upper  and  lower  bounds  exist  such 
that  k{  <  R(Y,T,Y,T)  <  k2  along  some  segment  of  the  ray,  then  the  period  of  affine 
parameter  between  consecutive  conjugate  points  satisfies  nj y[k^  <  T  <  nf y[k^  .  If  the 
sectional  curvature  is  zero  or  negative  everywhere  along  yF  {X) ,  then  neighboring  rays 
will  eventually  diverge.  (A  single  caustic  may  form  if  the  initial  region  is  concave  but 
otherwise  one  does  not  expect  the  periodic  development  of  focal  points  further 
downstream.)  In  situations  where  a  ray  occasionally  may  enter  convergence  and 
divergence  zones,  one  cannot  conclude  from  checking  the  Riemann  tensor  whether 
conjugate  points  will  occur  and  how  frequently. 
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Appendix  B 


NULL  GEODESICS  AND  ACOUSTIC  RAYS 


This  section  introduces  the  connection  between  the  null  geodesics  of  a  pseudo- 
Riemannian  manifold  and  the  bicharacteristics  of  the  equations  of  hydrodynamics  of 
ideal,  isentropic  fluids  for  the  case  of  subsonic  flow.  The  emergence  of  a  geodesic 
structure  from  the  study  of  hyperbolic  partial  differential  equations  is  a  hallmark  of  the 
method  of  characteristics.  The  same  geodesic  structure  emerges  in  the  study  of  acoustic 
rays,  a  connection  that  has  been  pointed  out  in  the  past  by  researchers  in  both  the 
acoustics  and  relativity  communities.  In  most  cases  one  appeals  to  the  high-frequency 
approximation  although  it  should  be  pointed  out  that  the  same  structure  emerges  in  the 
general  treatment  of  the  system  describing  the  propagation  of  discontinuities.  The 
majority  of  this  section  contains  little  new  information  and  is  written  for  completeness, 
rigor  and  convenience,  because  the  variety  of  different  approaches  has  not  been  cataloged 
in  one  place.  Every  attempt  at  completeness  is  made  while  maintaining  brevity. 

B.l  The  Bicharacteristics  of  the  Equations  of  Hydrodynamics 

The  equations  of  hydrodynamics  for  ideal,  isentropic  fluids  consists  of  the  continuity 
equation  and  Euler’s  equation,  respectively  given  by 


dtp  +  o  ■  Vp  +  pc2V-u  =  0  , 

(B.la) 

+  +  —  Vp  +  V®  =  0, 

(B.lb) 

P 


where  p  is  the  fluid  pressure,  v  the  fluid  velocity,  p  the  fluid  density,  and  where  the 
fluid  obeys  an  equation  of  state  p  =  p(p) ,  with  p'  =  c~2  defining  the  local  sound  speed. 
The  potential  ®  is  due  to  external  forces  such  as  gravity  acting  on  the  fluid.  Equations 
(B.la)  and  (B.lb)  are  a  system  of  first-order  quasi-linear  partial  differential  equations. 
Applying  the  method  of  characteristics  to  this  system  leads  to  the  characteristic  matrix 

Dcp  pc2i$<pf  ,  (B2) 

_p^(p  KiD(P  \ 

where  the  function  cp(x,t )  is  the  wavefront,  defined  as  a  surface  in  space-time  on  which 
the  initial  data  of  the  fields  are  specified. 
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The  characteristic  conditions  are  satisfied  by  setting  the  determinant  of  the  matrix 
equal  to  zero  leading  to 

Q  =  (D(p)%D(p)2 -c2V(p-V(p)=  0.  (B.3) 

Each  term  in  Eq.  (B.3)  defines  a  partial  Hamiltonian  for  the  system  and  may  be  treated 
separately.  Defining  Hl=(D<pY  and  H2  ={^D(p)2  -  c2V  cp  V  (p),  one  set  of 
characteristics  is  determined  by 

H2=gMVdM<pdv<p  =  0,  (B.4) 


where  the  contra  variant  metric  tensor  guv  has  been  introduced  with  g00  =  —c  2 , 
g0'  =  -pi ,  g‘J  =  Sj  -  ni/ij ,  with  ju  =  o/c  defining  the  local  Mach  number,  or  in  matrix 
form 


-1  -y, 

C2S  ;;  -0:0 


K~VJ 


J  'J 


(B.5a) 


The  corresponding  covariant  metric  tensor  given  by 


r-c2+o 2 


S„y  = 


v  -°J 


8 


Ji  J 


(B.5b) 


is  the  inverse  of  Eq.  (B.5a)  satisfying  g uagav  =  8V .  A  factor  of  pc 2  has  been  removed 

from  Eq.  (B.4)  for  simplicity  but  otherwise  having  no  effect  on  the  geometry  of  the 
curves  defined  by  it.  The  bicharacteristics,  which  are  identified  with  the  curves  along 
which  energy  is  propagated,  or  with  the  rays  of  the  system,  are  given  by  a  set  of 
trajectories  in  space-time.  Defining  an  appropriate  parameter  X  along  these  curves  and 
introducing  the  canonical  momentum  p p  =  d  (p  in  //2 ,  the  bicharacteristic  equations  are 


dxfl 

dX 


1  dH2 

2  dPf, 


gMVPv’ 


(B.6a) 


dpM  _  1  8H2 

dX~  2  dx » 


\paPpdnS 


a(3 


(B.6b) 


In  Eq.  (B.6b)  the  derivative  of  the  contravariant  metric  can  be  replaced  by  the  derivatives 
of  the  covariant  metric  by  differentiating  g flclgav  =  8V ,  thus  combining  Eqs.  (B.6a)  and 
(B.6b)  leads  to  the  following  equation  for  the  bicharacteristics: 
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which  is  precisely  the  equation  for  the  geodesics  of  a  differentiable  manifold  M  with 
metric  g  .  From  Eq.  (B.7)  and  the  characteristic  condition,  it  follows  that  the 

bicharacteristic  curves  are  equivalent  to  the  null  geodesics. 

B.2  Linear  Acoustics 

To  describe  acoustics  in  the  presence  of  background  fluid  motion,  all  fields  in  the 
problem  are  written  in  the  form  o  =  w+v ,  p  =  pQ  +  px,  and  p  =  p0  +  px,  where  w ,  p0 

and  p0  are  respectively  the  velocity,  pressure  and  density  of  the  fluid  medium,  and  v , 
px  and  px  are  contributions  due  to  acoustic  phenomena.  In  the  following  treatment  of 
the  problem  it  is  assumed  that  the  perturbations  are  weak,  leading  to  a  linear  theory  for 
the  propagation  of  sound.  No  restriction  is  placed  on  the  state  equation  for  the 

background  fluid  motion,  while  the  acoustic  fluctuations  are  assumed  to  obey  the 

adiabatic  condition,  px=c~2px.  The  background  fields  are  assumed  to  obey  the 

equations  of  hydrodynamics;  with  body  forces  such  as  gravity  lumped  together  in  the 

background  equations,  independently  of  the  acoustic  field.  Applying  these  assumptions 
to  Eqs.  (B.la)  and  (B.lb)  gives 

\d*  px  +  c2p0V  ■  v}+  [pxV  ■ w  +  c2v  ■  V p0  +  pxc2 D* c  2 1 =  0 ,  (B.8a) 

\p0D*v  +  Wpx }+  \c~2  pxD*w  +  PqV  ■  Vw}  =  0,  (B.8b) 

for  the  acoustic  perturbations,  where  D  =  8t  +  w-V ,  following  Thompson’s  notation 
and  the  adiabatic  assumption  has  been  used  to  eliminate  derivatives  of  px .  It  is  not 

assumed  that  the  background  quantities  or  their  variations  are  small.  Letting  the  four- 

component  vector  u[,]  denote  the  field  variables  with  u[,]  =vx  for  i  =  2,  3,  4,  and 

u [1]  =  px ,  Eqs.  (B.8a)  and  (B.8b)  may  be  written  in  matrix  form 

Lm[u]  +  M-u  =  0,  (B.9) 

where  Z(1)[---]  is  a  first-order  differential  operator  acting  on  u  while  M  is  an  ordinary 

matrix  acting  on  u.  In  both  cases  these  matrices  depend  on  the  background  fields  and  the 
acoustic  fields: 

L  Jb*  W!(vf  ]  m\c-1V-w  +  D'c-2  (Vpj 
Lv  K,P»D'\  [  C~2D‘  W  p0V®iv_' 


47 


One  can  see  by  inspection  that  the  first  tenns  in  Eqs.  (B.8a)  and  (B.8b)  resemble  Eqs. 
(B.la)  and  (B.lb),  leading  to  the  conclusion  that  the  characteristic  condition  for  acoustic 
perturbations  in  the  presence  of  background  fluid  motion  is  given  by  Eq.  (B.5a),  where 
the  terms  in  the  metric  contain  the  background  w  instead  of  v  . 

B.3  The  Eikonal  Approximation,  Acoustic  Rays 

Applying  D  to  Eq.  (B.la)  after  dividing  through  by  pc2 ,  taking  the  divergence  of  Eq. 
(B.lb),  and  taking  the  difference  of  the  resulting  equations,  leads  to 

D-^—Dp-V-Vp  +  A,  =0,  (B.lOa) 

pc  P 


while  taking  the  gradient  of  Eq.  (B.la)  and  applying!)  to  Eq.  (B.lb)  leads  to 

DpDu -Vpc2V  ■  o  +  A2=0,  (B.lOb) 


where  A l=-dkujdjok  and  A 2  =-9kpVuk .  Applying  the  same  procedure  to  Eqs. 

(B.8a)  and  (B.8b)  leads  to  a  similar  result  for  the  linear  acoustic  field  in  the  presence  of 
background  motion: 

\d*-^D*Pi  -V— VPl  l  +  At  =0,  (B.lla) 

[  Poc  Po  J 

{fi)*/70Z)*v  -  V/?0c2V  •  v}+  A2  =  0,  (B.llb) 


where 


Aj  =D 


f  \ 

Av-# 

P0c 


-V 


-A D*w  +£>*(v  -V\np0)+D* 

Poc  J 


A  D  ^ 

—  D*c~2 
\Po  J 


-5*(v-VwJ-5*vAw*> 


A2  =  D*{pov  ■  Vw)+  D*(pxD*w)~  v(pjV  •  w)-  v(c2v  •  Vp0 )-  8kpxVwk . 


The  equation  for  the  eikonal  in  the  high-frequency  limit  of  linear  acoustics,  like  that  for 

the  bicharacteristics,  comes  from  the  highest-order  derivatives.  The  terms  At  and  A2 
contain  only  the  field  variables  and  first-order  derivatives.  The  eikonal  approximation  is 
defined  by  taking  a— >0,  pl=7iae"p,a  and  v  =  aae"/,la,  with  na  =  7r0  +  aKx  +  •  •  • , 
da  =  <r0  +  ad  j  +  •  •  • ,  where  complex  fields  have  been  used  for  simplicity.  Inserting  these 


48 


expansions  for  the  field  variables  and  collecting  all  terms  of  order  a  2  leads  to  the 
following  eikonal  equation  in  matrix  form: 


~-H2 

or" 

0 

B 

(B.12) 


where  H2  is  defined  in  Eq.  (B.3)  and  B  =  c2V cp  0  V<p-  \^{p*  <p) Defining 
Q  =  —H2  ©  B  for  the  block-diagonal  matrix  in  Eq.  (B.12),  the  condition  for  the  existence 
of  a  nontrivial  solution  to  Eq.  (B.12),  detQ  =  ~H2  detf?  =  0 ,  leads  to  the  condition 

Q  =  {p*(pfi^D*(pf  - c2V (p-'S/ (f^j  =  0.  (B.13) 

By  inspection  of  the  terms  Aj  and  A2 ,  one  sees  that  they  may  be  set  to  zero  in  the  limit 
of  a  slowly  varying  environment  since  every  term  contains  derivatives  of  the  background 
fluid  variables.  This  leads  to  a  decoupling  of  the  equations: 

D*  — ^—D*px  -  V  •  —  V/?j  «  0 ,  (B.14a) 

Poc  Pq 

D* p0D*v  -  Vp0c2V  •  v  «  0  .  (B.14b) 


B.4  Time-Parameterized  Rays  in  Three  Dimensions 

The  surface  of  constant  time  (p{x,t)t=h  =  cp(x) ,  generated  by  taking  a  constant  time 

cross-section  of  the  wavefront,  is  referred  to  as  a  phase  surface  (Fig.  Bl).  (From  here  on, 
a  wavefront  is  a  surface  of  equal  time,  and  “in  phase”  is  taken  to  mean  equal  time  of 
flight.) 

In  general,  displacements  in  the  null  hypersurface  are  subject  to  the  constraint 
g00dt2  +  gydx'dx1  +  2 g0idtdx'  =  0 ,  which  leads  to  the  following  constraint  on  dx/dt : 

g^dx_  dx^  +  2g^dx_  =  _1  (B  15) 

kon  dt  dt  gim  dt 

For  the  specific  metric  given  in  Eq.  (B.5a),  Eq.  (B.15)  states  that  the  ray  velocity  relative 
to  the  moving  fluid  equals  the  local  sound  speed  {dx/dt  -  o)- (dx/dt  -  fi)  =  c2  (see  Fig. 
B2). 
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Figure  Bl.  Sample  of  a  complete  space-time  ray  structure.  (A)  Family  of  geodesics  labeled 
by  s,  where  the  s  parameterized  flow  lines  are  constructed  from  the  deviation  vector  acting 
on  a  fiducial  geodesic.  (B)  Each  flow  line  represents  a  curve  of  constant  affine  parameter.  A 
tangent  vector  and  the  deviation  at  a  single  point  are  illustrated.  (C)  A  curve  of  equal-time 
(wavefront)  constructed  from  the  intersection  of  a  constant  time  plane  with  the  geodesic 
flow.  Equal-parameter  curves  do  not  lie  in  the  constant  time  surface  in  general. 


V 


Figure  B2.  Relationship  between  the  fluid  velocity, 
ray  tangent,  and  wavefront  normal  in  Euclidean  space. 
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An  equation  for  dx/dt  is  derived  from  Eq.  (B.7)  by  writing  the  equations  for  time  and 
space  separately: 


d  x  k  dx‘  dxJ  k  dt  dt  k  dxl  dt 

- r - h  I  ij - hi  00 - h  21  (0 - =  0  . 

d/X  d/ 1  d/ 1  d/ 1  d/ 1  d/ 1  d/ 1 


(B.16a) 


d  t  o  dxl  dx'  o  dt  dt  0  d.xl  dt 

-  +  1  ij - hi  00 - h  21  ;0 - =  0  . 


d/J 


d/ 1  d/ 1 


d/ 1  d/ 1 


d/ 1  d/ 1 


(B.16b) 


Changing  the  parameter  along  the  ray  from  X  to  t  leads  to 


dx 


d  x  _  d  t/dX  dx'  _  ^k  dx'  dxJ  _  rk  _  ^k  _  _A 

- z h  — - — - h  1  ij - hi  00  +  21  /  0  -  —  U  . 

dt  (dt/dX)  dt  dt  dt  dt 


(B.17a) 


d2t/dX2 

(dt/dX)2 


0  dx'  dxJ  o  ^OT.o  dx' 

+ 1  ij  — : - : — 1-1  oo  +  Z1  w  — —  =  U  . 


dt  dt 


dt 


(B.17b) 


Using  Eq.  (B.  17b)  to  replace  the  second  term  in  Eq.  (B.  17a)  gives 


d2xk 
dt 2 


0  dx'  dx1  o  o  dx' 

1  ij - hi  00+21  i  0 - 

dt  dt  dt 


dxk 


dt 


k  dx  dxJ  k  k  dx1 

+  1  ij - 1-1  00+21  i  0 - =  U  . 

dt  dt  dt 


(B.18) 


Using  the  Christoffel  symbols,  Appendix  D,  and  some  algebra,  the  equation  for  the  3- 
dimensional  ray  paths  as  a  function  of  time  becomes 

At  =  2^ sij (A  - U )(u  -  Vj )  +  fe  ~  vk td-o  Inc  +  2r  -  Vine) 


2c 


+ 


8tuk-  cdkc  -  (okj  (xj  -  Vj )  ■ + 1  (Sv  -  mkj  )uj . 


(B.  19) 


Equation  (B.19)  is  converted  to  an  equation  for  the  wavefront  normal  through  the 
identification  cn  =  x-o  ,  and  can  be  used  to  determine  either  the  evolution  of  x  or  h  : 


+  +nklfil  =0, 

dt 


(B.20) 


where  Q  kl  =  1  sjki  (v  x  v  +  s  x  n .  +  2Vc  x  h)  . . 
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Along  with  a  complete  presentation  of  the  bicharacteristics  and  a  discussion  of  the  variety 
of  situations  described  by  null  geodesics,  the  main  results  of  this  appendix  are  Eqs.  (B.7), 
(B19),  and  (B.20).  In  particular,  note  that  the  geodesics  determined  by  Eq.  (B.7)  describe 
the  propagation  of  discontinuities  and  the  evolution  of  the  wavefront  in  the  most  general 
situations,  for  nonlinear  disturbances  propagating  in  a  random  environment.  The  catch- 
22  is  that  in  order  to  find  these  curves  in  the  nonlinear  case  one  needs  the  complete 
velocity  field,  o  =  v  +  w ,  which  requires  knowing  the  complete  solution  in  the  first  place. 
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Appendix  C 


CONSTRUCTION  OF  THE  INTERNAL  FRAME  FIELD 


C.l  Ray-Centered  Frame  in  Space-Time 

To  calculate  the  geometric  spreading  of  an  acoustic  beam  from  the  deviation  equation 
the  Riemann  tensor  must  be  projected  into  a  ray-centered  frame  field.  The  frame  being 
defined  formally  on  the  null  hypersurface  requires  a  basis  of  four  mutually  orthogonal 
vectors  in  space-time.  Following  the  notation  in  Hawking  et  al?2,  the  ray  tangent,  T M ,  is 
chosen  to  serve  as  one  of  the  coordinate  directions.  A  second  null  vector,  If ,  satisfying 
the  condition  LfTvg  =  -1 ,  is  chosen  as  the  second  basis  vector.  To  complete  the  local 

geodesic  coordinate  basis,  two  space-like  directions,  and  ef ,  satisfying  the  conditions 
T"e/a  =  0 ,  LaeIa  =  0 ,  and  e“eJa  =  8 u  for  I,  J=  1,2,  are  introduced. 

The  two  null  vectors  define  a  time-like  direction  Et  =  c_1(l  fi)  that  is  orthogonal  to 
the  two-dimensional  space-like  hypersurface  (e^  at  all  points  in  space-time.  After 
defining  this  pseudo-orthonormal  basis  at  an  initial  point  on  yF  (A) ,  the  basis  is  parallel 
transported  along  yF(A)  to  define  a  new  basis  at  each  point  by  solving  the  parallel 
transport  equation: 

+  Ta/tvT,t(eI)v  =0-  (C.l) 

dA 

The  null  vector  If  is  automatically  parallel  transported  along  the  null  geodesic  T . 

To  ensure  that  neighboring  rays  are  emitted  from  the  source  at  the  same  “time”  as  the 
fiduciary  ray,  the  initial  space-like  basis  is  chosen  to  be  purely  spatial  in  the  coordinate 

frame,  (e, )"  =  (()  e, ) ,  and  the  initial  deviation  T0"  =  (o  Y0 )  is  chosen  to  lie  along  one 
of  the  initial  basis  vectors. 

Equation  (C.l)  formally  preserves  the  constraints  placed  on  the  initial  conditions. 
Beginning  with  16  quantities,  4  vectors  with  4  components  each,  the  conditions  impose 
11  constraints  leaving  5  independent  quantities.  Four  of  these  are  the  components  of  T  “ 
determined  by  the  geodesic  equation,  leaving  only  one  remaining  quantity  to  calculate. 
The  vector  LM  may  be  determined  algebraically  in  terms  of  the  environmental 
parameters.  The  remaining  degree  of  freedom  may  be  identified  with  a  two-dimensional 
rotation  of  the  space-like  basis  {eI }  in  the  {eI )  hyperplane  expressed  as  R9ueJ() . 
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C.2  Ray-Centered  Frame  in  Cartesian  Space,  Wavefront  Normal  Basis 

In  this  section  an  alternate  choice  of  internal  basis  is  presented.  This  representation, 
tenned  auxiliary  basis,  is  better  suited  for  visualization  and  interpretation  of  the 
deviation,  as  it  always  remains  tangent  to  the  surface  of  constant  phase  in  space.  From 
the  conditions  TMe]  g^v  =  0  and  efe^g  =  1 ,  one  may  verify  that 

h -{e,  -  e°jt)=  0 ,  ||e7 -fe°||  =  1 ,  (C.2) 

with  respect  to  the  ordinary  dot  product  in  three-dimensional  space.  Hence  the  set  of 
three-dimensional  vectors 


~k  _  ~k  -o  dx 

ei  —  ei  ei 


dt 


(C-3) 


is  tangent  to  the  wavefront  everywhere  along  the  ray  and  the  basis  (h,  er )  remains 
orthonormal  as  it  is  propagated  along  the  ray.  The  vectors  e,  are  referred  to  here  as  an 
auxiliary  basis.  An  equal-time  deviation  is  constructed  by  projecting  the  components  Y, 
onto  the  auxiliary  basis.  The  deviation  projected  into  the  different  basis  will  point  in 
different  directions  in  four  dimensions  but  the  magnitude  is  independent  of  the  choice  of 

=  jFrV 


basis, 


Y 


2+Y{ 


The  sectional  curvature  needed  for  the  spreading  equation  may  be  rewritten  as 
4 Kjj  =Rftvap A/'vAja/i,  where  A/1'  =Tt‘e\-Tve>ll.  The  mixed  components  of  this 
tensor, 


Aok  =  To~k_Tk~o  =  To^k  _^dxk /dt))  =  T%k , 
are  clearly  proportional  to  the  auxiliary  basis.  The  pure  space  components, 


(C.4a) 


Al/C  rri  l  k,  rri  k  *  k  rji  l '  '  k  rri  k.  r^J  l 

i  =T  eI  -T  e1  =T  eI  -T  e,  , 


(C.4b) 


are  equal  to  the  components  of  the  cross  product  T  x  e,  in  three  dimensions.  The  tensor 
Ay'u,/  is  parallel  transported  along  y(X) , 

— A/v  +Tf,apTaAIpv  +TvapTaK^  =  0.  (C.5) 

dA 

Equation  (C.5)  tracks  six  fields,  counting  indices  in  four  dimensions.  It  is  clear  from  the 
definition  of  A/1  that  there  are  only  three  degrees  of  freedom.  By  separating  Eq.  (C.5) 
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into  individual  components,  the  following  transport  equation  for  the  auxiliary  basis  e/  in 
Ei  is  derived: 


—  e,  +  —  (v  x  o  +  s  x  n .  +  2  Vc  x  n ) ,  el  =  0 , 
dt  2  1  y 


(C.6) 


with  sk  =  njSjk  .  Note  that  an  overall  factor  ofdt/dA  has  been  divided  out.  In  principle, 

the  parallel  transport  equation  has  a  solution  in  terms  of  a  path-ordered  integral  involving 
the  SO(3)  generator  evaluated  along  the  ray.  In  practice,  these  are  not  known  in  closed 
form.  However,  assuming  that  the  ray  equation  has  been  solved  either  analytically  or 
numerically  the  nonnal  n  is  known,  reducing  the  problem  to  a  two-dimensional  rotation 
about  the  wavefront  normal  followed  by  an  application  two  Euler  matrices  to  orient  the 
basis  along  the  ray  path. 

Choosing  the  initial  basis  to  match  a  global  Cartesian-coordinate  basis 
(n(0)  e,  (0)  e2(0 )}  =  (?'  j  k") ,  the  rotation  matrix  taking  /7(0)  — >  /7  may  be  found 

explicitly.  The  rotation  of  (e7  \  about  h  is  denoted  R(o) .  The  action  of  this  sequence 
on  the  initial  basis  vectors,  interpreted  as  an  active  transfonnation  on  the  initial  values 
(n(0)  e,  (0)) ,  is  denoted  U  =  R((p)R(6)R{a )  .  The  rotation  R(a)  has  a  block-diagonal 
form  R(a)  =  1  ©  r  ,  where  r  is  a  two-dimensional  rotation  about  the  x-axis  of  the  global 
coordinate  system.  Inserting  this  into  Eq.  (C.6)  leads  to  the  following  differential 
equation  for  r: 


dr  „ 
- b  (OCT  r 

dt 


0, 


where 


<7,. 


' 0 
v1 


-n 


( 0  = 


nCl 


yz 


2  2 

n;  +  n; 


The  solution  to  Eq.  (C.7)  is 

^cosa  -sin  o' 
r  = 

sin  a  cos  a 


(C.7) 


(C.8) 


where  a  =  J  a>\  dt ,  with  the  integrand  evaluated  along  the  specific  ray  path  in  space.  The 

dependence  of  a>  on  specific  Cartesian  indices  comes  from  the  choice  of  lining  up  the 
initial  basis  with  the  global  coordinate  system. 

Equation  (C.8)  reduces  the  problem  of  tracking  eight  degrees  of  freedom  plus 
constraints  to  one  of  evaluating  a  single  integral.  This  reduction  is  not  without  cost.  It  is 
clear  that  the  denominator  will  vanish  whenever  a  ray  has  a  horizontal  turning  point. 
While  this  may  not  happen  often  in  modeling  underwater  systems  where  attention  is 
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focused  on  parabolic  rays  trapped  in  sound  ducts,  it  will  happen  in  more  general 
situations.  The  danger  of  this  divergence  occurring  indicates  that  it  may  be  more 
practical  to  track  the  components  of  the  larger  dynamical  system,  using  the  constraints  to 
regulate  error  along  the  way. 

The  ability  to  reduce  Eq.  (C.l)  to  Eq.  (C.7)  depends  on  the  identification  of  A^r  with 
the  wavefront  normal  basis.  In  general,  this  trick  may  always  be  employed  to  reduce  the 
number  of  quantities  needed  to  calculate  the  sectional  curvature  and  is  invaluable  in  2- 
dimensional  ray  systems,  3-dimensional  space-time  systems,  as  it  allows  one  to  generate 
e  by  rotating  h  by  90  .  Without  this  insight  one  would  be  compelled  to  either  track  all 
three  components  of  e  by  Eq.  (C.l)  or  evaluate  Eq.  (C.8). 
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Appendix  D 


RIEMANN  TENSOR  AND  CHRISTOFFEL  SYMBOLS  FOR  THE  ACOUSTIC 

METRIC 


For  the  acoustic  metric  defined  by  the  characteristic  condition,  Eq.  (B.5), 


1 

r-i 

r-c2+^2 

c2 

1-°/ 

°25ji  -°jvi) 

9  ^  //V 

l  "yy 

^  J 

direct  substitution  into  Eqs.  (A.4)  and  (A.  10)  lead  to  the  following: 


and 


r°oo  =  d  v  Inc  + 

2  c2  ,J 


r°!0=ain c--%s,  , 

'  2 c2  9 


r'oo  =  Otd_0  In  c  -  dp,  +  cG,c  ~  7;^%  ~  (c28ki  -  op,  )sjk  , 


F'yo  =  ofi j  In c  +  \Q)j~L^r sjk  . 


rC  -  n  r°  •  -  k  9 

1  J  -  uk 1  V  -  2c2  ’ 


D  _  ^  C  _  C  C  ) 

' nimj  A  2  \  mn^  ij  ^  mi nj  ) 

4  c 


(D.la) 

(D.lb) 

(D.lc) 

(D.ld) 

(D.le) 


(D.2a) 


Rm,  =  \(8,S„  ^c-S^,  lnc)  +  ^(^„ -VJ  ■  (D-2b) 

Rmj  =-^SlSl)-pXajSm  +  8j‘»l,)+c8l8JC  +  ^S:j8,  lnC  ~LjXLS"Sf 

+  y  (V<  lnc  +  S*3;  Inc  -  ^.0.  Inc)  -  ±(s„S*  +  ®WS„  +  ®*sj  . 

(D.2c) 
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Typically  the  Riemann  tensor  contains  second  derivatives  of  the  metric  components  as 
well  as  first  derivatives  squared.  A  significant  feature  of  this  system  is  the  lack  of 
second-order  time  derivatives  of  the  environmental  parameters  and  the  coupling  of  dtc  to 

the  fluid  velocity.  For  a  stationary  medium  the  focusing  is  determined  solely  by  the 
spatial  derivatives  of  the  sound  speed. 
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Appendix  E 


FERMAT’S  PRINCIPLE  AND  SNELL’S  LAW 


E.l  Fermat’s  Principle 


Fermat’s  principle  follows  immediately  from  the  line  element  by  solving  for  dt: 


dt  = 


■J(v  -  dr)2  +  (c7  -  dr  -  v  ■  dr 


2  2 

c  -v 


(E.l) 


where  the  +  root  is  taken  for  subsonic  flow.  For  supersonic  flow,  dr  ■  o  >  0  and  the 
travel  times  are  detennined  by  the  -  root. 

E.2  Metric  Symmetry  and  Snell’s  Law 

An  isometry  (symmetry  of  the  metric)  is  indicated  by  a  cyclic  coordinate30.  The 
presence  of  a  cyclic  coordinate  in  the  metric,  labeled  xc ,  leads  to  a  conservation  law  for 
the  corresponding  momentum  or  conjugate  variable, 

dPc  dg^T*  d&g^e  q 

dX  dX  dX 


The  vector  called  a  Killing  vector,  points  in  the  direction  of  the  xc  coordinate 

curves.  (The  statement  “...presence  of  a  cyclic  coordinate”  xc  is  equivalent  to  the 
absence  of  that  coordinate  from  the  metric  altogether.)  When  there  is  enough  symmetry, 
the  second-order  geodesic  equations  are  reducible  to  a  set  of  first-order  equations. 

Specializing  to  a  time-independent  environment,  the  Killing  vector  =  (l  6)  leads 
to  the  conservation  law 


/  2  dt  _  dr 
-C  -u~  —  -O-—  =  Kn 


’  dX 


dX 


(E.3) 


Traditionally,  this  is  associated  with  the  energy  of  the  particle,  and  for  the  metric 
signature  in  Eq.  (B.5)  k0  is  negative.  (From  here  on,  k{]  >  0  and  an  overall  minus  sign  is 

dropped  from  Eq.  (E.3).)  Consider  the  problem  of  a  layered  medium  where  c  and  0 
depend  on  only  one  Cartesian  coordinate,  z.  There  are  two  more  conserved  currents, 
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which  may  be  thought  of  as  components  of  the  ray’s  translational  momentum. 
Killing  vectors  i  and  j  lead  to 


dr  dt  _ 

- oT  =  k  , 

dX  dA  T 


The 


(E.4) 


where  k  =  Kxi  +  K2j  and  uT=  oj  +  ovj .  Imposing  the  null  constraint  on  the  components 
of  the  tangent  vector  leads  to 

-c2i2 +k-k  +  {z-u,i\  =  0  (E.5) 

and 


c2i  =  k0-ut  ■  k -uz{z -uzt).  (E.6) 

Equations  (E.3)  and  (E.4)  may  be  used  to  derive  a  generalized  version  of  Snell’s  law 
originally  presented  by  Kornhauser14  for  the  special  case  vT  =  vi  .  By  definition  Eq. 
(E.3)  states  that  k  =  ichT  which  in  turn  states  that  the  projection  of  the  wavefront  normal 
in  the  x-y  plane  points  in  a  fixed  direction.  Equation  (E.3)  may  be  rewritten  to  give 
k{)  =  f  ■  ch  =  ic(c  +  v  ■  n) ,  relating  the  energy  to  the  angle  between  the  wavefront  normal 

and  the  spatial  component  of  the  ray.  When  u  =  oT  =  vi ,  Eq.  (E.6)  becomes 
tc0  =  ic(c  +  on x )  .  Taking  the  momentum-to-energy  ratio  in  this  case  gives  Snell’s  law  for 
the  wavefront  normal: 


cos  6 

- =  a  , 

c  +  v  cos# 


(E.7) 


with  a  =  ATj / at0  and  nx  =  cos 6* .  The  ray  angle  is  related  to  9  by  ccosO  +  o  =  rcosy/  . 
The  magnitude  of  the  ray  velocity  follows  from  the  null  constraint 


r  =  ucosi//  + 


-  v2  cos°  y/  . 


(E.8) 


Insertion  of  Eq.  (E.8)  into  Eq.  (E.7)  gives  Kornhauser’s  result  for  Snell’s  law: 


1 

a  =  — 

c  1 


£  sin 2  y/  +  cosy/yjl  -  s 2  sin2  yr 
s 2  sin2  y/  +  cosy/^l- s 2  sin2  y/ 


(E.9) 


The  paths  of  the  rays  propagating  in  a  layered  media  with  =  0  are 
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(E.10) 


dt  _  [k()  -ut  -k) 
dA  c 2 


dr  o  /  _  _ 

—  =  k  +  — r l/q,  -  vT  •  K 
dX  c2  V  0  r  ' 


dz  _  -J(/r0  -  vT  ■  k)2  -  K  ■  KC2 

dA  c 


(E.ll) 


(E.12) 


Taking  appropriate  ratios  of  Eqs.  (E.10)  and  (Ell)  leads  immediately  to  the  ray-path 
integrals  encountered  in  underwater  acoustics. 

The  procedure  may  be  generalized  to  any  set  of  coordinates.  A  derivation  of  two- 
dimensional  Snell’s  law  in  polar  coordinates  is  now  presented.  The  acoustic  line  element 
in  polar  coordinates  is 

-(c2  -v2)dt2  -2urdrdt  -2verd6dt  +  dr2  +r2d8 2  =  0.  (E.13) 


Specializing  to  a  rotationally  symmetric,  time-independent  environment  with  o  =  vO , 
leads  to  the  following  set  of  equations  for  the  ray: 


2  2*2  2  \  2 

err  =  -a~  +  (r  -  aoy , 

(E.14) 

c2r20  =  ac2  +  v{r  -  ao )  , 

(E.15) 

a 

8 

1 

5s. 

II 

<N 

(E.16) 

Following  all  the  same  steps  leads  to  the  same  result  for  Snell’s  law,  Eq.  (E.9),  with  0 
interpreted  as  the  instantaneous  angle  between  the  wavefront  normal  and  9 . 
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